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AESTRACT 


The  cnject  of  this  thesis  has  osen  twofold.  The 
first  object  was  to  develop  FORTRAN  versions  of 
several  existing  APL  programs  which  were  designed  tc 
analyze  univariate  data.  In  particular  the  programs 
were  designed  to  test  for  exponentiality  and  normality 
of  the  data  and,  by  sectioning  or  jacknifing,  obtain 
estimates  of  sampling  variances  of  sancie  moments. 
The  second  object  of  the  thesis  was  tc  use  these 
programs  in  a  simulation  study  of  first-ccne 
first-served  queues  in  which  the  service  times  and  the 
inter-arrival  times  are  exponentially  distributed  but 
dependent.  The  dependence  is  introduced  by  using  the 
mixed  moving  average  autoregresive  structure 
(  EAR23A(p,q)  )  for  exponential  sequences  introduced  by 
Lewis  and  cc-workers.  Four  models  cf  correlated 
queues  are  introduced,  giving  autoccrrelated  and 
cross-corieiate d  service  and  arrival  times  in  various 
degrees.  The  simulation  study  gives  a  quantitative 
idea  cf  the  effect  cf  correlation  en  the  mean  waiting 
time  and  the  distribution  of  the  waiting  time. 
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I.   INTRODUCTION 


The  object  of  this  thesis  has  been  twofold. 

The  first  object  was  to  develop  FCETHAN  versions  of 
several  existing  APL  programs  which  were  designed  to  analyze 
univariate  data,  in  particular  to  test  fcr  exponentiality 
and  normality  of  the  data  and,  by  sectioning  or  jackrifing 
obtain  estimates  of  variances  of  sample  moments. 

The  second  object  was  to  use  these  programs  in  a 
simulation  study  of  first-come  first-served  queues  in  which 
the  service  times  and  the  inter-arrival  times  are 
exponentially  distributed  but  dependent.  The  dependence  is 
introduced  by  using  the  mixed  moving  average  autoregr essive 
structure  (  EAflMA(p,g)  )  for  exponential  sequences 
introduced  by  Gaver  and  Lewis  (1977),  Lawrance  and  Lewis 
(1977),  Jacobs  and  Lewis  (1977)  and  Lawrence  and  Lewis 
(1978).  Four  schemes  are  introduced  which  give 
autocorr elated  and  cross-correlated  service  and  arrival 
times  in  various  degrees. 

Since  no  analytical  properties  of  these  queues  can  be 
derived,  their  properties  are  studied  ty  simulation. 
Because  the  EABMA  queues  are  neither  regenerative  nor 
aarkovian  there  are  several  problems  in  the  simulation  study 
of  knowing  when  one  can  assume  that  a  steady  state  has  teen 
reached  in  the  simulation  and  of  obtaining  confidence 
intervals  for  estimated  parameters.  This  problem  is  handled 
by  generating  a  number  of  independent  sample  paths  and 
comparing  statistics  obtained  from  the  replications  at 
several  points  along  the  sample   paths.    Only   the   waiting 
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time   process  (W  )  is  studied  in  this  thesis,   not  the  state 

N 
space  process  N (t)  . 


Unfortunately  the  simulation  study  of  the  EAEM5-type 
queue  was  hampered  by  the  size  and  speed  of  the 
IBM  System  360/67  system  at  the  Naval  Postgraduate  School. 
Consequently  the  results  were  not  as  extensive  as  had  teen 
hoped  for.  Fairly  detailed  results  are  given  for  the  case 
where  the  service  time  process  is  autoregressive ;  one  case 
where  the  traffic  intensity  is  0.99  and  the  correlation 
parameter  0.98,  was  intensively  investigated.  This  gives 
seme  idea  of  the  length  of  the  transient  in  the  queue  and  of 
the  inflation  of  the  mean  waiting  time  because  of  the 
correlaticn  in  the  service  process.  Some  investigation  of 
the  case  where  the  service  times  and  inter-arrival  times  are 
cross-correlated  are  given. 

In  section  II  we  describe  the  basic  histograi  and 
statistics  computation  package  'HISIGS/F'  and  apply  it  to 
t*o  sets  of  telephone  error  data.  In  section  III  the 
plotting  method  and  the  plotting  Subroutines  'NGFHEL', 
' BXPLT'  are  described  and  we  test  them  with  various 
generated  data.  In  section  IV  we  describe  the  assessment  of 
variability  Subroutines  'S2CTN',  'JACK'  and  we  apply  them  to 
telephone  data  set  1.  In  section  V  the  Subroutine  'LIST*  is 
described  and  we  apply  it  to  two  sets  of  telephone  error 
data.  In  section  VI  we  describe  and  analyze  the  EAEM5(p,q) 
model  as  well  the  program  which  we  developed  in  order  to 
simulate  it.  Computer  program  listings  are  provided  after 
section  VI. 
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II.   SUBROUTINE  HISTGS/HISIE S 


A.   DESCRIPTION 


Here  the  Subtoutine  'HISTGS/HISTFS  •  is  presented. 
•HISTGS/HISIf S '  is  used  for  estimating  the  probability 
density  function  from  a  given  set  of  data  and  at  the  same 
time  computing  some  basic  statistics.  Easically  this 
Subroutine  is  the  library  FCRTRAN-Subroutine  HISTF/G,  which 
was  created  at  N.P.S.  by  D.  H.  Robinson.  It  has  teen 
modified  ty  adding  the  new  entry  point  'LIMITS1. 

The  new  entry  point  may  be  used  for  wild  data  01  for 
data  having  a  mixed  distribution.  Therefore  it  is  a  useful 
tool  for  that  user  involved  with  this  kind  cf  data.  Simply, 
the  user  has  to  define  those  values  in  the  range  of  data 
(sectioning  the  data  into  several  disjoint  sections)  that 
he/she  relieves  are  useful  cut  points.  For  each  section  of 
data  the  user  may  have  the  probability  density  function  as 
well  as  the  basic  statistics  (  as  HISTF/G  dees  ) .  Also  by 
sectioning  the  data  over  its  range,  the  user  may  have  a 
histogram  cf  the  random  variable  when  it  is  conditioned  to 
te  between  given  limits. 

The  number  of  sections  that  the  entry  pcint  'LIMITS*  can 
accept  is  as  many  as  5C. 

A  complete  description  of  how  'HISTGS/HISTFS'  operates 
is  given  in  the  subroutine.  On  the  other  hand,  a  summary  of 
the   subroutine   is   given   by   typing   on   the  terminal  the 
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command  DESCEIEE  HISTGS  under  CMS  (Cambridge  Monitor  System, 
flef.  4).  Iyping  'DESCRIBE  HISTGS'  the  following  response 
will  be  given  on  the  terminal: 

SUBRCQTINE    HISTGS/HISTES 

' HISTGS/EISTFS ■  gives  you  a  histogram  of  a  set  of  data 
along  with  the  estimate  of  a  set  of  basic  statistics.  Also 
it  gives  you  the  following  options: 

1.  Displaying  a  smoothed  empirical  probability  density 
function  ever  the  histogram, 

2.  Sectioning  the  data  into  several  disjoint  sections 
(no  more  than  50)  and  then  having  a  histogram  and  the 
estimate  of  tasic  statistics  for  each  section, 

3.  Scaling   the  histogram, 

4.  Displaying  just  the  histogram  (with  or  without  the 
density  function)  and  no  statistics. 

The  estimated  basic  statistics  are:  Mean,  Median, 
Irimean,  Midrange,  Geometric  Mean,  Harmonic  Mean,  Variance, 
Standard  Deviation,  Coefficient  cf  Variation,  Mean 
Deviation,  Range,  Midspread,  Third  and  Forth  Central  Moment, 
SXewness,  Kurtosis,  Beta  1 ,  Beta2,  Maximum,  Minim  u  it,  and 
Cuantiles . 

CALLING   SEQUENCES 


CALL 

HISIGS 

(  X,  N, 

NBAR  ) 

CALL 

HISIfS 

(  X,  N, 

NEAR  ) 

CALL 

LIMITS 

(  X,  N, 

XLIM, 

CALL 

NCSIAS 

CALL 

SIAIS 

CALL 

FIX5    ( 

SCALE  ) 

CALL 

NCEIXS 

CALL 

PECMAX 

(  PSC  ) 

CALL 

NCPRMX 
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ARGUMENTS 
X      Array  of  data  values 
N      Number  of  data  values 
NEAR    Number  of  bars  in  the  histogram 
X1IM    Array  of  cut  points 
LIMS    Number  of  cut  points 

SCALE   Vector  of  two  values  to  scale  the  data 
PSC    A  real  variable  (between  0.  and  1.)   defining 
the  maximum  value  of  the  probability  axis. 

If  the  entry  point  NOSTAS  is  called  before  callinc  any 
of  the  entry  points  HISTFS,  HISTGS,  or  LIMITS  then  no 
statistics  are  printed,  otherwise  the  statistics  are  printed 
ty  default. 

By  calling  the  entry  point  FIXS  the  histogram  scale  may 
be  fixed  and  remains  set  unless  it  is  reset  ty  another  call 
to  FIXS  or  allowed  to  float  (the  default)  by  calling  NCFIXS. 
The  reason  for  fixing  the  scales  is  so  that  when  comparing 
mere  than  one  batcn  of  data,  we  can  get  ccuparably  scaled 
histograms. 

By  calling  PBOMAX,  the  maximum  value  of  the  probatility 
axis  (Y)  can  be  set  to  PSC,  which  will  remain  in  effect 
until  NOPEMX  is  called  (the  default) . 

More  information  is  given  in  tne  sunroutine. 


E.   USING  'HISTGS/HISTFS'  WITH  TELEPHONE  DATA  1  AND 
TELEPHONE  DATA  2 


Here  '  HISIGS/HISTFS'  has  been  used  to  plot   the   set   of 
Telephone   Data   1   and   Telephone   Data  2.   The  algorithmic 
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procedure  fcr  both  sets  of  data  is  the  same.   That   is,   the 

entry   point   'HISTFS'   is   called   first  and  then  the  entry 

point  'LIMITS'.   Pour  cut  points  have  been  used  (1,  2,   141, 

86000)   sectioning   the  data  into  three  disjcint  sections  as 

fellows: 

Section  1  :   Consists   of   data   points   X    ,  such  that 

<i> 

1  <  X     <  2. 

<i> 

Section  2  :   Includes   the   data  pcints   X    ,  such  that 

<i> 

2  <  X     <  141. 

Section  3  :   Includes   the   data  points   X    ,  such  that 

<i> 

142  <  X     <  86000. 
<i> 

No   print-out   for   section   1   is   given  because  of  the 

constant  value  of  the  data  pcints. 


Figures  1  through  6  give  the  histogram  and  the  tasic 
statistics  fcr  both  sets  of  data  (  For  the  entire  set  and 
fcr  sections  2  and  3  ) . 

We  observe  that  both  of  data  sets  are  so  unruly  that  the 
use  of  the  entry  point  •LIMITS'  cf  'HISTGS'  is  suggested. 

Observing  Figures  1  and  4  we  may  conclude  that  both  data 
sets  appear  to  be  the  same,  having  seme  geometric  or 
exponential  distribution.  But  as  we  will  see  later  this  is 
net  true.  Thus  the  histogram,  in  that  case,  which  is 
obtained  by  using  just  the  entry  point  'HISTFS'  is  not 
helpful  and  leads  us  to  a  wrong  conclusion. 

• LIMIIS'  is  applied  now  to  give  us  a  useful  answer  : 

Observing   figures   2  and  5,  obtained  by  using  the  entry 

point  'LIMITS',  we  may  conclude  that  the  data   sets   dc   not 

have   the   same   distribution  (  at  least  in  that  interval  ). 

In  addition  data  set  2  is  far  from  having  a  geometric  or   an 
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exponential  distribution,  as  it  can  be  seen  ty  the  histogram 

of  the  entire  set  of   data.    This   is  a   conclusion   which 

cannot   he   cttained   ty   without   use  of   the   entry  point 
•LIMITS' . 

Observing  also  Figures  3  and  6  ,  where  the  histogram  of 
the  3rd  section  from  each  data  set  is  presented,  we  may 
conclude  that  this  section  of  both  data  sets  may  have  an 
exponential  distribution,  but  not  the  same  one. 

The  most  important  facet  of  the  data  which  appears  is 
the  modes  at  anout  124.  The  mode  is  much  mere  apparert  in 
Telepnone  Data  2.  In  fact  closer  investigation  reveals  that 
this  mode  was  due  to  leakage  of  dial  pulses;  this  leakage 
caused  bit  errors  at  multiples  of  124  bits. 

The  above  is  a  quick  and  informal  analysis  of  the 
Telephone  Data  Sets  obtained  by  using  the  entry  point 
•LIMITS'  of  • HISTGS/HISTFS' ,  and  the  conclusion  is  that  both 
of  data  sets  are  so  unruly  that  any  formal  characterization 
of  their  distribution  is  not  so  easily  obtained.  In 
addition  it  is  apparent  that  a  mixed  distribution  governs 
both  sets  of  data. 

Having  quickly  analyzed  both  of  the  telephone  data  sets, 
the  use  of  entry  point  'LIMITS'  cf  the  Subroutine 
•aiSTGS/HISlPS'  is  obvious. 
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TELEPHONE    DATA    1 


THE    ENTIRE    SET    OF    EATA 
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Figure    2    -         TELEPHONE    DATA    1.       2nd    SECTION. 
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III.   PLOTTING  SUBROUTINES 


k.       INTRODUCTION 


Graphical  methods  of  assessing  the  validity  cf  a 
probability  model,  and  of  estimating  some  fcasic  parameters 
(especially  location  and  scale  parameters)  have  been  widely 
applied  over  many  years.  A  variety  of  different 
prescriptions  have  been  advanced  for  the  plotting  positions, 
that  is,  for  the  set  of  values  at  vihich  the  ordered 
observations  in  a  sample  should  be  plotted. 

Plotting  the  ordered  data  is  an  informal  and  quick 
visual  methcd  for  getting  an  impression  of  hew  the 
distribution   of   data   looks.    Depending   en  what  is  known 


about   the   distribution   of 
probanility  plotting  methods. 


the   data 


there 


are 


many 


The  Plotting  Subroutines  presented   here,   'NORMPL'   and 
•EXPLT1,   use   the   methcd   cf  plotting  the  crder  statistics 

(  X  )  versus  their  expected  value  (  III  ]  )  called 
•SCORES1.  This  method  does  not  require  knowledge  cf  the 
Cumulative  Distribution  Function  F  (x)  (except  for  its 
assumed  continuity).  Therefore  it  may  he  very  attractive 
for  those  whe  knew  nothing  about  the  data  tc  be  analyzed.  In 
addition  the  indication  of  the  graph  may  lead  the  analvst  to 
further  useful  results. 
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Using  the  Plotting  Subroutines  'NORMPL'  and  lEX£LXa  cne 
may  analyze  the  resulting  graphs  obtained  tc  get  an  idea  of 
the  distribution  of  the  data. 


B.   DESCBIPTICN  OF  THE  PLOTTING  METHOD 


The  idea  behind  this  methcd,  as  abcve  stated,  is  tc  plot 

X     versus  E[  X    ]. 
<i>  Ci> 


It   is   used   where   data   X  ,  X  ,    .  .  .  ,  X     arise 

1    2        n 


as 


independent   observations   of  a  continuous  Bandom  Variable  X 

with  a  distribution  function  which  is  believed  to  have   seme 
particular     form     F (x)  .     Then     the    ordered    data 


x   ,  x 
C1>    C2> 


..,  X     are  plotted  against  Efx    1. 

CQ>  Ci) 


If   the   model  for  F(x)  used  to  get  E[ X    ]  is  true  then 
the  resulting  plot  is  a  linear  (regression)  relationship  and 

should   plot,   within   limits   of   sample   variations,   as  a 

straight  line.  Barnett  [2]  gives  the  value  cf  E[ X    ]   as   a 

C  i> 

linear  relation  of  the  plotting  position  x    .   That  is: 

Ci) 

EfX    ]  =  m  +  sx        , 
<i>  Ci) 

here  m,  and  s,  are  the  location  and  scale  parameters.  The 
Mean  and  Variance  of  X  are  net  necessary.  hs  we  see  the 
above  relation  uniquely  defines  the  plotting  position  as  : 
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EC  X     -  m] 

.JiL 

<i)        s 


X     -  m 
<i> 
=  EC ] 


Therefore,  if  the  model  is  correct  the  linearity  will 
probably  be  evident.  In  addition  if  the  resulting  plot  not 
being  linear  leads  us  to  reject  the  model,  we  may  get  an 
indication  cf  the  type  cf  departure  present  and  thus  an 
indicaticn  cf  what  action  tc  take  or  what  alternative  mcdel 
to  use.  It  should  be  noted  that  the  linearity  of  the  graph 
suggests  informally  that  the  model  is  true,  but  net  tc 
accept  the  mcdel.  Additional  formal  tests  must  be  done  in 
crder  to  accept  the  model. 


C.   SUBROUTINE 


NCRMEI 


1 .   Description 

Subroutine  'NCRMPL',  is  used  tc  test  a  sample  of 
independent  observations  of  a  continuous  Random  Variatle  X 
having  a  distribution  function  which  is  believed  to  have  the 
Normal  form. 

'NOREEL'  sorts  the  data  into  increasing  crcer  to 
create  a  set  cf  Order  Statistics.  Then  the  Crder  Statistics 
are      plotted     versus      either      Normal      Scores 


(if   n  <  50,   n  =  number  of  data-points)  or  inverse 


the 


-  1 


Standard   Normal  Cumulative  Distribution  Function   (F   (y) ) , 
(if   n  >  50)  . 
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Evaluation   of   Normal   Scores  is  based  on  computing 

th 
the  Expected  Value  of  the  i   Order  Statistic   (  Sf X     )  . 


This  evaluation  follows: 


Let 


th 
g (x    )   be   the   density   function   cf   the   i    Order 

Statistic  of  a  sample   X  ,  X  ,  %..,  X    from   a   population 

1    2        n 

having   the  Standard  Normal  Distribution  F.   Then  g  (x    )  is 

<  i) 

given    by    : 


g  u       ) 


n»  i-1  n-i 

F      (X  )  (1-F      (X  ))  f      (X  ) 

(1-1)'       (n-i)»        X       <i>  X       <i>  X       <i> 


n» 


(i-1)  I      (n-i)! 


i-1      n-i 

p         Q         z      , 


where: 


Z    = 


f     <*,-J 
X       <l> 


-x  /2 

<i> 


(2  r\) 


V2 


*ttj 


P    = 


F     <x        ) 
X       <i> 


7 


i«o(2  n) 


1  -u   /2 

e  du 

V2 


Now 


Q    =    1-E     . 


fOO 


E[  X         ]    =      x    g  (x)    dx 
<i>         J 

-co 


(1) 


CC 


n"  i-1      n-i 

= x    Z    P  Q  dx 

(i-1):    (n-i);  J 
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The  value  of  expression  (1)  has  been  evaluated  and 
is  given  in  [19],  This  table  gives  the  values  of  Normal 
Scores  frcm  1  to  n/2  (for  n<50)  .  The  values  from  (n/2)+1  to 
n  are  evaluated  using  the  symmetric  expression  : 


Also 


B[X  ]    =    -    E[X  ]       . 

<i>  <n-i+1> 


E[X  ]    =    0,  n   =    2k+1       , 

(k+1) 


•NOEaEL'  uses  the  above  Normal  sccres  if  n  <  50.  For 
n  >  50  • NQfiaEI*  computes  its  own  Normal  Sccres  using  the 
program  function  'INVNRM',  which  gives  the  value  x  cf  the 
inverse  Standard  Ncrmal  distribution  (x  =  F  (y) , 
0  <  y  <  1) .  Justification  of  this  is  the  fact  that  if  the 
data     are     really     normally     distributed     then 

E[  X    ]    F(  i/ (n+ 1)  ).    Therefore   •NORMPL'  can  accomodate 

<i  > 
any  sample  size. 

In  addition  to  the  plot  given  by  •NORMEL',  the 
Wilk-Shapiro  test  for  normality  can  be  evaluated  by  the  same 
Subroutine.  A  derailed  description  of  that  test  is  given  in 
[22].  However  a  summary  of  that  test  is  presented  here  in 
order  to  help  us  to  describe  how  the  program  works.  The 
Wilk-Shapiro  test  is  based  on  the  statistic: 


«  =  —2   . 

S 


where : 


b  = 


k 
1=1  n-i  + 1 


(X 


(n-i  +  1) 


k=r/2,    n  e\en 
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b  =  a  (X    -X    )  *  a    (X      -X    )  +  .  .  . 
n   <n>   <1)      n-1   <n-P   <2> 


+  a    (X      -X    )  , 
k+2   <k  +  2>   <k> 


k=  (n-1)  /2,   n  odd 


2    n         2 

s  =  f:  (x  -  i) 

i=  1   i 


lh€  coefficients  (a  )  are  defined  by: 

i 


n 
a   =  SI  m  u   /C  , 

i    j  =  1  j  ij 

where : 


i  -  1,2,...,n 


m   =  E[ X    ], 
j       <j> 

u   =  CcvfX    ,X    ], 
ij       <i>   <j> 


(That  is  the  'NORMAL  SCCRE51) 
i#  j  =  1  ,2,  .  .  .  ,n 


2    ■  -1  -1 

C  =  j  V   V   m 


where 


id  =  ( m  ,  m  , .  .  .  m  ) 
1   2     n 


V  =  (u    ),    That  is  the  nXn  Covariance  Matrix. 


Ther€   are   some   approximations  associated  with  the 
evaluation  of  these  coefficients.   However  fcr   n   <  5C   the 

(a  )   are   tabulated   in   the   above-mentioned  paper  and  the 
i 

program  uses  them  by  FCBTRAN  DATA  initialization  statenents. 
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For   n>50  the  (a  )  are  computed  by  the  program  using 

i 
the  following  approximation  method: 

a   =  2m  /C,        i  =  1,2,...,n-1   , 

i     i 

where  C  is  given  by  the  following  Least-Sguare  equation  as  a 
function  of  n: 

2 

C   =  -2.722  +  4.C83n 


Also,  for  n=1 


2      P(n/2  +  .5) 
a   = 

2    f(n/2  +  1) 


Note  that    a  =  a 

i    n-i+1 


Seme  properties  for  the  W  Statistic  can  be  given 
here,  which  have  neen  taken  from  the  above  mentioned   paper: 

(a)  .   W   is  scale  and  origin  invariant. 

(b) .  8  has  a  distribution  which  depends  only  en  the 
sample  siza  n,  for  samples  from  a  Normal  Distribution. 

2 

(c) .   m   is  statistically  independent  of  S   and  of  1      for 

samples  from  a  Normal  Distribution. 

r       2r     2r 
(d)  .   E[W  ]  =  E[  b   ]/E[S   ]   for  any  r. 

(e) .   The   maximum  value  of  I  is  1  . 

2 
(f)  .   The   linimum  value  of  H  is  na  /(n-1)  . 

(g) .   The   half  and  first  moments  of  U  are  given  by: 

2 
r  i/2      B  Rn/2  -.5) 
E[  H    ]  =  

V2 
C  p<n/2)  2 

2   2 
B  (B  +1) 
E[W]  = 

2 
C   (n-1) 
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where: 


2    •  -1 

B  =  n  7   m 


(h)  .   For   n=3  the  density  of  W  is  given  as: 

3        -1/2  -V2 

t     (w)  =  ~z-    (1-w)  w    ,  .75  <  n  <  1 

H        n 


2.   Procjiam S  tructure 

■NORBPL*  is  a  FOBTR AN-callable  Subroutine  with  each 
call  returning  (optionally)  either  a  plot  fci  a  given  set  of 
data  or  a  value  for  the  W  Statistic  or  both. 

The  program  is  divided  into  three  parts: 
The  first  part  is  the  control  prcgraa  of  'NOFMEL'. 
In  that  part  the  'Normal  Scores'  for  a  =  2(1)25  and 
n  =  26(2)50  are  stored  along  with  the  Subroutine  by  CATA 
statements;  on  the  ether  hand,  the  Normal  Sccres  for  r.  >  50 
are  computed  by  calling  the  program's  function  'IN7NKM'.  In 
this  part  the  user  specifies  just  a  plot  cr  just  the  value 
of  W  or  tcth.  The  Library-Subroutine  'PXSGEI'  is  used  by 
the  program. 

The  second  part  of  the  program  consists  of  the 
Subroutine  'PLOT'  which  accepts  any  data  set  to  be  tested 
for  Normality.  'PLOT'  itself  scales  the  data  according  to 
its  range  and  plots  the  scaled  data  along  110  equal  spaced 
positions  of  the  X-axis.  No  plot  is  given  if  the  data  has  a 
constant  value. 


The   third   part   of   the   program   is   the    'WILK1 

subroutine.     In   this   Subroutine   the   coefficients   (a  ) 

i 

(for  n  <  50)  of  the  H  Statistic   are   stored   using   FCBTEAN 


36 


DATA   initialization   statements.   For  n>50  the  coefficients 
(a  )  are  computed  by  the  Subroutine  itself.   Ihis  Subroutine 

also   calls  the  'INVNRM'  Function  to  get  the  'Normal  Scores' 

to  be  used  fcr  the  computation  cf  the  (a  )  's  vihen  n  >  50. 

i 

A  complete  description  of  how  'NOEHEL'  operates  is 
given  in  the  Subroutine.  However  a  summary  is  obtained  by 
typing  on  the  terminal  DESCRIBE  NOH21PL .  When  the  user  types 
the  command  DESCRIBE  NCRHPL  under  the  CMS  environment  the 
following  response  is  printed  on  the  terminal: 

SUBROUTINE    NCRHPL 

'NOBaEL'  takes  a  set  of  data,  sorts  it  into 
ascending  order  and  uses  the  created  Order  Statistics  for: 

1.  Plotting    X     versus   either   'Normal   Scores'   {if 

(i) 

q  <  50,   n  =  Sample  size)   or  Inverse   cf   the   unit   Normal 

Cumulative  Distribution  (that  is: 

-1 
X     versus  Normal  scores  or  F   (  i/(n+1)  )   ,) 
(i) 

to  see  if  there  is  a  linear  fit. 

2.  Computing    the   value   of   the    H    test    Statistic 
(Milk-Shapiro  tast  fcr  normality) . 

It  is  called  by: 
CALL    HCBMEI   (  X,  SCCRES,  N,  K  ) 

ARGUMENTS 
X      Is  the  array  containing  the  data 
SCORES  Is  a  work  array  of  dimension  N 
N      Is  the  number  of  data  values 
K      User's  option 

For:    K  =  1     A  plot  only  , 
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K  =  2    The  W-value  is  given  , 

K  =  3     Both  of  the  above  are  given. 

Mere  information  is  given  in  the  subroutine. 

3 •   Interpreting  The  Output 

If  data  are  really  Normally  distributed  then  a 
straight  line  is  expected  from  the  output.  This  is  an 
indication  that  gives  the  user  a  first  feeling  of  the 
distribution  of  data.  But  if  the  plot  is  not  linear  then 
the  user  has  to  reject  the  normality  assumption.  In 
addition  the  shape  of  the  plot  may  lead  the  user  to  make  an 
alternative  decision  of  what  model  he  may  use.  For   example: 

(a) .  Saape  of  the  form  of  Figure  7a  suggests  that 
data  are  skewed  to  one  side  (specifically  positive  values). 
In  that  case  an  Exponential  model  could  he  an  appropriate 
one.   The  histogram  of  the  data  is  given  in  Figure  7b. 

(b) .  Shape  of  the  form  of  Figure  3a  can  be 
interpreted  as  data  not  as  dispersed  as  the  Normal 
Distribution.  This  set  of  data  may  have  a  Unifcru  cr  a 
Triangular  Distribution,  and  could,  for  instance,  arise  if 
measuring  tolerances  or  components  hand-picked  to  lie  within 
tolerances.   The  associated  histogram  is  given  in  Figure  8b. 

(c) .  Shape  of  the  form  of  Figure  9a  indicates  much 
more  dispersed  data  than  Normal,  but  the  data  is  symmetric, 
protably  coming  from  a  Double  Exponential  or  Cauchy 
Distribution.  The  associated  histogram  is  given  in  Figure 
9b. 

It  should  be  noted  here  that  a  straight  line  does 
not   always   arise   from  Normal  data.  Figure  10a  for  example 
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gives  an  indication  that  data  are  Normal  (Straight  lice). 
But  these  data  have  been  generated  by  computer  from  a 
Symmetric  Triangular  Distribution  in  the  range  (0  to  2).  On 
the  other  hand  figure  11a  plots  a  data  set  cf  size  50  again 
from  a  Triangular  (0  to  2)  distribution.  Eut  in  this  case 
we  observe  that  the  graph  is  not  a  straight  line.  These 
results  demonstrate  the  informality  of  that  Probatiliry 
Elotting  Method  and  in  these  cases  the  user  has  to  continue 
testing  the  data  by  other  formal  tests.  Nevertheless,  since 
this  method  is  primarily  an  initial  screening  device,  ease 
cf  application  is  important.  Figure  11b  gives  the 
associated  histogram. 
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Figure    7a    - 


GENERATED    DATA    FROM    AN    EXPONENTIAL 
DISTRIBUTION 
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Figure    7b    -  HISTOGRAM    OF    THE    DATA    Or    FIGaRE    7a 
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Figure  8a  - 


GENERATED  DATA  FROM  A  UNIFORM  (0,1) 
DISTRIBUTION 
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Figure    8b   -         HISTOGRAM    OF    THE    DATA    CF    FIGURE    8a 
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Figure  9a  - 


GENERATED  DATA  FROM  A  CAUCEY  DISTRIBUTION 
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Figure    9b    -  HISTOGRAM    OF    TH2    DATA    CB    FIGURE    9a 
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Figura  10a  -    GENERATED  DATA  PROM  A  TRIANGULAR  (0,1) 
DISTRIBUTION.  SAMPLE  SIZE  200 
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Figure    10b    -         HISTOGRAM    CF    THE    DATA    CJ    FIGURE    10a 
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Figure  11a  -    GENERATED  DATA  FROM  A  TRIANGULAR  (C,1) 
DISTRIBUTION.  SAMPLE  SIZE  50 
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Figure     11b    -  HISTOGRAM    OF    THE    DATA    Of    FIGURE    11a 
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4«   Using  'NORflPL1   with   Generated   Data   from   Various 
Distributions  . 


The  Subroutine  'HOfifllL*  has  been  used  here  tc  test 
for  Nomality  some  data  sets  which  were  generated  ty  the 
computer  frca  various  distributions.  The  plots  of  these 
data  sets  illustrate  the  advantages  and  disadvantages  cf  the 
plotting  Subroutine  •NCHMPL*. 

Observing  Figures  7a-b,  8a-b,  9a-b  and  12a-b  we  may 
note  the  following: 

(a) .  Figure  12a  shows  a  plot  of  a  sample  of  5C  data 
points  from  a  Normal  Distribution  (  N  (1000,1)  ).  As  it  was 
expected  a  very  straight  line  is  fitted.  In  addition  to 
that  plot,  the  Subroutine  'NORMPL'  uses  the  option  to 
evaluate  the  ■ HIIK-SHAEIRO1  test  and  the  W-value  is  confuted 
and  printed  en  the  same  figure.  Comparing  this  value  with 
the  percentage  points  cf  the  H-test  which  are  given  in  [22], 
for  sample  size  n=50  we  see  that  we  may  accept  the  norsality 
assumption  with  a  significanse  level  =  .C2  .  Figure  12b 
gives  the  associated  histogram. 

(t) .  Figure  7a  plots  a  sample  cf  200  data  points 
cf  the  Unit  Exponential  distribution.  Here  the  nonlinear 
graph  is  obvious  (  as  is  expected  )  and  suggests  tc 
reject  the  normality  assumption  of  the  data  without  doing 
any  further  formal  test. 

(c) .   Figure  8a   gives  a  plot  cf  a   sample   cf   size 
m 
200   from   Uniform   (  C,  1  )     generated  data.  The  shape  of 

that  plot  indicates  data  not  as   dispersed   as   the   nerval, 

therefore   the   user  is  told  to  reject  normality,  getting  at 

the  same  time   a   suggestion   that   data   may   be   Uniformly 
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distributed   (  because   of   that   particular   shape   of   the 
graph  ) . 

(d) .  Figure  9a  suggests  strongly  departure  frcm  a 
Normality  assumption.  Indeed  this  set  of  data  have  teen 
generated  from  a  Cauchy  Distribution  and  its  shape  really 
indicates  a  symmetric  distribution. 
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Figure  12a  - 


GENERATED  DATA  FROM  A  NORMAL  (100C,1) 
DISTRIBUTION 
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Figure     12b    -  HISTOGRAM    OF    THE    DATA    OP    FIGURE    12a 
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5.   UsiUS  'NORMPL*  With  Cost  Overruns  Data 

'NOEHPL'  was  used  with  the  ccst  overruns  data  [6], 
for  the  year  1950  to  see  if  there  is  a  linear  fit.  Locking 
at  Figure  13  we  see  that  data  may  be  Normally  distributed 
since  the  plct  is  very  close  to  linear.  Of  course  scaie  data 
points  deviate  from  the  straight  line,  but  these  points  may 
be  considered  as  outliers. 
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Figure    13    - 


COST    OVERRUNS     LATA 
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D.   SUBROUTINE  EXPIT 


1 .   Description 

Subroutine   'EXPLT'   is   indended   tc   be   usee   for 

testing  a  saiple  of  independent  observations  cf  a  continuous 

random   variable   X   with   a   distribution  function  which  is 

telieved  to  fce  Exponential,  by  plotting  each   of   the   Order 

Statistics  X     versus  its  Expected  Value,  and  by  estiaating 
U) 

the  parameters  ^  and  Y  of  the   observations   (  y  ,   and  \) 
introduced  later) . 

This  Subroutine  sorts  the  data  into  increasing 
order,  obtaining  the  Order  Statistics.  The  Order  Statistics 
are  then  plotted  versus  Exponential  Scores.  These  scores  are 
evaluated  by  the  Subroutine.   The  Exponential  Scores  are  the 

expected  value  cf  the  Order  Statistics   (  E[ X    ]  ) .     Ihey 

(i) 
are  more  easily  derived  and  computed  than  the  Normal  Scores. 

That  is: 

Let  X  be  a  random  variable  Exponentially  distributed 
with  parameter  >=  1 
Then  it  is  keewn: 


-^xx 

f  (x)  =  ^€   ,      a,  x  >  o  , 


F  (x)  =  1-e     ,        V  x  >  0 

A 
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th 
Ihe  Density  cf  the  i   Order  Statistic  is  given  by: 

n-1  i-1  n-i 

f     (x    )  =  n(      )  F  (x    )     (1-F  (x    ))     f  (x    ) 
X      <i>       i-1    X   <i>         X   <i>       X   <i> 
<i> 

st 
Thus  for  the  1    Order  Statistic  we  have: 

-  *x         -;\x 

<*>    n-1      tn 
£     (x(i))  =  n(1-(1-e      )  )    >e 

(1) 


-n>x 

<  1 ) 


=  n.\e 


Therefore  X     is   exponentially   distributed   with 
(1) 


parameter  ">  =  n  >  . 
And  for  *\  =    1  w<=  have: 


1      1 

E[X    ]  = ■  =  — 

(1)      n\  n 


It  can  be  proved  (see  Feller  [7])  that  the  n   random 

variables   X    ,  X    -X    ,  X    -X    ,  ...,  X    -X       from 
(1)    (2)   (1)    (3)   (2)         (n)   (n-1) 

an  Exponential  random  variable  X  with  densitj  function: 

f  (x)  =  >e     ,         x  >  0    , 
X 

are  independent  and   the   density   of   the   random   variable 

X  =  X    -  X      is  given  by: 
(*)    (k-1) 


-  ^(n-k+1)  y 
f  (y)  =  (n-Jc  +  1)  *e  ,      y  >  0 


That   is,   Y   is  an  Exponential  randcm  variable  with 
parameter   (c-k+1)  y  . 
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Therefore    E[Y]      is    given   by: 


1                           1 
E[x         -X  ]    = = , 

(k)        (k-1)  ^(n-k  +  1)  n-k+1 


for        =1 


ihis    fact    is    used    to    derive    E[X         ],  i    =    2,3,.. .,n    . 

(i) 

lhat   is: 

E[X         )     =    E[      (X         -X  )     +     (X  -X  )♦...+ 

(i)  (i)       (i-1)  (i-1)        (i-2) 

(X  -X  )      +    X  1 

(2)        (1)  (1) 


Sf  X  ]    +    E[X  -X  ]    +    ...+ 

(1)  (2)        (1) 

E[X  -X  ]    +    E[  X  -X  ] 

(i-1)         (i-2)  (i)        (i-1) 


111  1 

n  n-2+1  n-3+1  n-i+1 


=  X  —2— 

k=1       n-k+1 


i    =    2  , .  .  .  n 


•EXPLT'      derives      the    values   of      E[ X         ]        using   the 

(i) 

anove  formula  and  then  it  uses   its   Subroutine   ' EPLC1*   tc 
plot  the  Crdsr  Statistics  versus  their  expected  value. 


In  addition  to  the  plot,  which  is   an   informal  and 

guick    test    for    exponentiality   the   estimates   of  the 

statistics  V  and  /f^  provide  another  informal  test  since  the 

values    of  i     and  V   have    a    constant    value    for  any 
Exponential  distribution: 
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£y    definitioc     ^jLand    Y     are    given   by 


Ihus 


a 


(J4 


cr3 


13  2  2         3 

—    £[X      -    3X  y    *   3Xy    -  y.  ] 


(1) 


13  2  2  3 

--     (Z[X     ]    -    3j*E[X     ]    +    3y  E[X]    -y  ) 


Now 


3 
E[X     ] 


08    3 
x      f    (x)     dx 

~     3  -*i 


■  I  *  y 


ye  dx 


/•"4-1         4       -;\x 

re*)         x        ^    ^       d* 


rco 


The  above  integral  is  that  of  the  Gamma  density  function 
with  parameter  OL  =  4 ,  and  since  its  value  equals  1  we 
have: 


3      3i 
E[X  ]  =  — ■ 


.A 
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2 

Computing      E[X    ]      we    have: 

2  100    2  -;u 

E[  X    ]   =  x       ^  e  dx 

rp)    r  x       ^   e      dX 

2 
"    "*2 

Also  B[X]  »   ~   -   ^ 

Substituting    in    (1)    we    have: 

3 

Ji«  ec — ±—  ] 

16                    12                         111 
=    ..        ( 3 +    3 } 

<?3        j3  *     y2  y2      1        23 

3  3 

=    >      (    2/ J,     )     =   2       . 

Therefore    the    value   of      V         equals      2      for      any      Exponential 
distribution . 

Also     x        is    computing    as   follows: 

4  3  2    2  3  4  4 

=    E[  X      -    4y»X       +    6  yX      -4y<X+y]/C~      "3 

4  3  2         2  3  4  4 

(2)  =     {E[X     ]    -    4aE[X    ]    +    6y  E[  X     ]    -    4U  E[X]    +    U    }     /<J-3. 
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Eut 


4  (°°    4  -}X 

E£X    ]    =      \    x      ^e  dx 


r(5) 

4 


oo       5-1       5      -ax 

x  >     e  dx 


T(5) 


=    24   /    > 


Thus 


1  24  16  12  11 


(J1*  }4  7        >° 

4  4 

=     1      (9/>)"3=      6 


^       V 


V       X 


1 


-    3 


Therefore        the        value      of      ^      equals      6       for      any 
Exponential    cistr ibuticn. 

This      fact  is    used    by   the    Subroutine    'EXPLT*    to    give 

the    user    one    more  informal      test      for      exf cnentiality .        It 

estimates           the  values           of   ^    and      ^L       and        if      L  t    2 

and/or     \x  #    6    then  the   data      may      not      have      an      Exponential 
distribution. 

Estimating      the      values      of  ^L  and    Y   ,      'EXPLT'       uses 

3  4 

unbiased      estimator      for     & ,      2[  (X-^)     ]      and    E[  (X-y)     ]    using 

the    formulas: 

n  2  1/2 

a  =    [X   (X.-X)    /    (n-1)    ) 
1-1      i 


EC  <X-f )      ]    = 


n  .3 

.^    (x.-   X) 

1=1      1 

(n-1)     (n-2) 


EC  (X-)<)      ] 


n  _    4 

(n(n-2)+3)     (X    (X.-X)    ) 
1=1       l 


(n-1)     (n-2)      (n-3) 


3  C  (n-1)     (2n-3) 
n     (n-2)      (n-3) 
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2.  11*1111    Structure 

The  Subroutine  'EXFIT'  is  also  a  ECRTRAN  callable 
subroutine  with  each  call  producing  a  plot  fcr  the  given  set 
of  data  and  at  the  same  time  calculating  and  printir.g  the 
estimated  values  of  the  parameters  1  and  V  of  the  data. 

Basically  the  program  is  divided  intc  two  parts: 

The   first   part   is   involved   with  all  required 

computations   calculating   the   Exponential  Scores  and  the 

various   statistics   needed    to    compute  the   V  and  }f 
parameters. 

The  second  part  is  the  Subroutine  'EELOT'.  'EEIOT' 
scales  tne  data  points  according  to  their  range  and  plots 
the  scaled  data  along  110  egual-spaced  positions  of  the 
X-axis.  No  plot  is  given  if  the  data  have  constant  value  or 
if  there  exists  a  data  point  less  than  zero. 

A  complete  description  of  hew  ' EXEIT'  operates  is 
given  by  the  Subroutine  and  a  summary  description  is  given 
on  the  terminal  by  typing  the  command  DESCSIEE  2XPLT  under 
the  OAS  environment.  When  the  user  types  DESCRIBE  EXELT 
the  following  response  is  printed  on  the  terminal: 


SUBROUTINE    EXPLT 

'EXEIT'  takes  a  set  of  data,  sorts  it  into 
increasing  order  and  plots  the  created  order  statistics 
versus  exponential  scores  (expected  values  of  order 
statistics)  computed  by  : 
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*[*,.,]  "  y(   i/(n-k+1)  ),    i  =  1,2,.. .,n   . 
(1)     *=1 

Alsc,   it   computes   the   estimates   of  H         and  ^ 
parameters . 


It  is  called  by  : 
CALL    EXELT   (  X,  SCOEES,  N  ) 

ARGUMENTS 
X      Is  the  array  containing  the  data 
SCOEES  Is  a  work  array  of  dimension  N 
N      Is  the  number  of  data  values 

Mere  information  is  given  in  the  subroutine. 

3 •   Interpreting  T he Output 

If  data  have  an  Exponential  distribution  then  the 
plot  will  tend  to  a  straight  line.  But  it  should  be  noted 
that  the  inverse  is  not  generally  true.  The  linearity  of  the 
graph  is  an  indication  only  and  gives  the  user  a  first 
feeling  of  the  distribution  of  data.  If  the  plot  is  nor 
linear,  however,  this  suggests  that  they  are  not 
Exponential.  In  addition  the  shape  of  the  plot  may  lead  the 
user  to  make  an  alternative  model  selection. 

Cn  the  other  hand  the  values  of  the 
statistics  jj(  and  jL  may  be  used  to  test  informally  the 
distribution  of  data.  These  values  are  an  indication  only 
for  the  test  for  expcnentiality.  If  the  estimate  of  ^  has 
the  value  about  2  and  the  estimate  of  ^ has  the  value  about 
6  then  it  suggests  that  data  may  be  exponentially 
distributed.  But  a  departure  of  these  values  suggests  that 
data  do  net  have  an  exponential  distribution. 
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**  •   Using   '  EXPLT1   with   Data   Generated   from   V ar ious 
Eistributions 


The  Subroutine  'EXPLT1  has  been  used  here,  plotting 
data  generated  by  the  computer  from  various  distributions, 
in  order  to  get  a  visual  sense  of  the  behavior  of  the   plot. 

Figure  14a  gives  the  plot  of  10C  Uniform  (0,1) 
variates.  enviously  ncnlinearity  governs  the  plot  (as  we 
should  expect)  and  thus  it  suggests  departure  from  the 
Exponential  distribution.  Besides  the  nonlinearity  of  the 
plot  the  estimates  of  ^and  ^  (0.1  and  -1.15  respectively) 
are  far  away  from  the  corresponding  values  (2  and  6)  of 
exponentially  distributed  data.  Figure  14b  gives  the 
associated  histogram. 

Figure    15a    plots    variates  generated   from   a 

Triangular  Symmetric  distribution  (0,2).  The   nonlinearity 

and   the   values   of  the  estimates  Y  and  jj   suggest  that  one 
reject  the  exponential  assumption. 

2 
Figure  16   shows   the   plot  of  X   generated  variates 

with  10  degrees  of  freedom.   For  the  same  reasons   as   above 
the  exponential  model  is  rejected. 

Figure  17a  plots  Normally  (N  (1000,1))  distributed 
generated  data  and  the  nonlinearity  of  the  plot  is  apparent. 
Eesides  the  nonlinearity,  the  estimates  of  ^L  and  ^ 
(.09  and  -0.63)  suggest  not  only  to  reject  the 
Exponentiaiity  of  data  but  in  addition  lead  us  to  test  the 
data  for  Normality.  Figure  17b  gives  the  associated 
histogram. 
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Figure  18  plots  exponentially  generated  variates 
with  parameter  ^ =  1.  As  we  see  the  graph  is  linear  anc  the 
user  gets  the  idea  that  the  data  possibly  are  exponential. 
Eesides  the  linearity  here,  the  values  or  the  estimates  of 
^  and  ^z  (2.06  and  5.4)  are  close  enough  to  their  true 
values  to  indicate  an  Exponential  model. 
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Figure  14a  - 


A  SAMPLE  CF  SIZE  100  FRCM  A  UNIFORM <0,1) 
DISTRIBUTION 
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Figure    14b    -         HISTOGRAM    OF    THE    DATA    Or    FIGURE    14a 
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Figure    15a    -  A    SAMPLE    OF    SIZE    100    FROM    A    TRIANGCLAR 

SYMMETRIC     (0,2)     DISTRIBUTION 
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HISTOGRAM    OF    THE    DATA    Of    FIGURE    15a 
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Figure    16a    - 


(10    d.f.)     DATA 
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Figure     16b    -  HISTOGRAM    OF    THE    DATA    CB    FIGURE     16a 
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Figure  17a  - 


A  SAMPLE  OF  SIZE  100  FROM  A  NORMAL  (1000,1) 
DISTRIBUTION 
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Figure    17b   - 


HISTOGRAM    CF    THE    DATA    Of    FIGURE     17a 
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Figure  18a  - 


A  SAMPLE  CF  SIZE  50  FROM  AN  EXPONENTIAL 
DISTRIBUTION 
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Figure     18b    -  HISTOGRAM    CF    THE    DATA    CP    FIGURE     16a 
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IV.   SUBROUTINE  LIST 


A.   DESCBIfTICN 


Another  Subroutine  presented  in  this  thesis  is  the 
Subroutine  fLIST'.  Its  purpose  is  to  list  a  given  set  of 
data  in  ascending  order  and,  taking  advantage  of  like 
occurrences  in  the  data,  to  print  the  ordered  data  in  a 
compressed  form.  This  feature  becomes  highly  useful  when 
listing  a  large  number  of  data  points  that  contain  many 
repeated  values.  It  is  also  a  tool  for  finding 
multiplicities  in  supposedly  continuous  data,  and  a 
probability  function  estimating  routine  for  data  which  is 
known  to  te  discrete. 

A  complete  description  cf  how  'LIST'  operates  is  given 
in  the  Subroutine.  However  a  summary  is  given  by  typing 
DESCRIBE  LIST.  When  the  user  types  the  command 
DESCRIBE  LIST  under  the  CMS  environment  the  following 
response  will  be  given  on  the  terminal: 


SUBROUTINE   LIST 

•LIST'  sorts  a  set  of  data   into   increasing   order   and 

gives  a  5-cclumn  print-out  as  follows: 

1st   column:   Serial   number  cf  the  first   occurrence   of 

this  value  in  the  ordered  List 

2nd   column:   Value  of  ordered  data-value 

3rd   column:   Frequency   cf  occurrence  cf  the  value 
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4th   column:   Percent   for  the  value 

5th   cclunn:   Graphical    representation  cf  the  frequency 
for  each  value. 


It  is  called  by: 

CALI    LISI  (  X,  N  ) 

where: 

X    Is  the  array  of  data 

N    Is  the  numner  of  data-values. 

Mere  information  is  given  in  the  subroutine 


E.   INT1BBHE1ING  THE  OUTPUT 


The  print-cut  is  a  visual  representation  cf  the  data  and 
of  each  data  value  frequency.  For  data  feints  having  the 
same  value  •IIST1  will  print  this  value  once  with  the  number 
cf  occurrences.  A  5-column  output  will  be  printed  and  its 
interpretation  is: 

1.  First    column   gives   the  serial  number  of  the  first 
occurrence  of  this  value  in  the  ordered  List 

2.  Second   column  gives  the  value  cf  ordered  data-value 

3.  Third   column  gives  the  frequency   cf   occurrence   of 
the  value 

4.  Forth   column  gives  the  percent  for  the  value 

5.  Fifth    column   is   a  graphical  representation  cf  the 
frequency  for  each  value. 
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^xamplej. 

Let  1,  1 ,    5,  3,  5,  1,  2,  2,  5,  6  be  a   given   data   set. 


Then 


the 


ordered 


data    will   be:    X   =X    =X    =  1, 

(1)   (2)    (3) 


X    =X    =  2,  X    =  3, 
<<*)   (5)       (6) 


X    =X    =X      5f 
(7)   (8)   (9) 


X     =6 
(10) 


Then 


the  print-out  will  have  the  following  form: 


~JKIJL_NUMBEi ORDERED  DATA FREQUENCY PERCENT PRO E^_ GRAPH 

.3 


1 

4 

6 

7 

10 


1 

2 
3 

5 
6 


3 
2 

1 
3 

1 


.  1 

.3 
.  1 


* 


If  there  are  no   data-ccints   having  a   common  value   then 

•LIST'   gives   this   indication   and  prints  ccly  the  ordered 

data.    This    happens    when    data  have    a  continuous 
distribution. 


USING  'LIST'  WITH  TELEPHONE  DATA  1  AND  TELEPHONE  EAIA  2 


Subroutine  'LIST'  was  used  with  Telephone  Data  1 
(Figure  1S)  and  Telephone  Data  2  (Figure  20)  and  a  brief 
analysis  of  the  output  follows. 

Looking  at  Figure  19  and  Figure  20  we  may  get  seme 
information  cf  each  data  set.  It  can  be  seen  that  both  of 
the  data  sets  contain  a  large  number  of  multiple  occurrences 
cf  the  data  value  one  and  data  value  two.  As  we  can  see 
reading  the  fourth  column  the  occurrence  cf  ones  is  1S9S  for 
Telephone  Data  1  and  24%  for  Telephone  Data  2. 
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Comparing  also  Figures  19a-c,  20a-d  we  can  see  that 
multiple  occurrences  happen  in  the  range  1  to  2UC  for 
Telephone  Data  1  and  in  the  range  1  to  132  for  Telephone 
Lata  2. 

Furthermore,  guick  visual  information  concerning  the 
rangs  where  *e  have  multiple  occurrences  can  be  obtained 
from  the  Probability  Graph.  Thus  we  can  see  that  for 
Telephone  Data  2  there  is  a  region  from  113  to  132  where  the 
multiple  occurrences  of  values  is  larger  than  in  neighboring 
regions.  Therefore  'LIST*  gives  the  user  a  useful  tool  for 
analysis  of  data. 


79 


SSRIAL  SUH3EH 


1 

129 
183 
211 
233 
250 
261 
271 
293 
297 
306 
316 
327 
333 
339 
345 
353 
361 
366 
378 
379 
384 
389 
392 
399 
402 
405 
407 
4  10 
4  15 

421 

425 
429 
431 

4J5 

4  38 

440 
442 
44  1 
445 
447 
448 
449 
453 
456 
457 
459 

4o0 
4t>  1 
4o3 

464 

465 

4b8 
470 
472 
474 
475 
477 
480 
481 
4d3 
484 

487 
468 
4S9 
490 
492 
49J 
496 
497 
498 
499 
501 
502 
504 
505 
506 


08DE8ED  DATA 


7BEQUENCIES 


0. 

10000000E 

01 

128 

0. 

2C0000002 

01 

54 

0. 

30000O00E 

01 

28 

0. 

UC0OOOOOE 

01 

22 

0. 

50000000E 

01 

17 

c. 

6COOO0OOE 

01 

11 

0. 

7C00O0O0R 

01 

10 

c. 

800O0000E 

01 

12 

0. 

90000000E 

01 

1U 

0. 

1C00OO00E 

02 

9 

0. 

1  100OOO0E 

02 

10 

0. 

120000002 

02 

1  1 

0. 

1 300O00OE 

02 

6 

0. 

140000JOE 

02 

6 

0. 

150U0OC0E 

02 

6 

0. 

loCuOOOOF. 

02 

8 

0. 

170010002 

02 

8 

0. 

18000OC0E 

02 

5 

0. 

19000000E 

02 

12 

c. 

2C000O00E 

02 

1 

0. 

210000002 

02 

5 

c. 

22000000S 

02 

5 

0. 

2300UOOOE 

02 

3 

0. 

24000000E 

02 

7 

0. 

25J0JG00E 

0  2 

3 

c. 

.'tiGUOOOOE 

02 

3 

0. 

27030000E 

02 

2 

u. 

2H00OOOGE 

02 

3 

y  • 

29000UOOS 

02 

5 

0. 

3C300000E 

02 

6 

0. 

3  100JOGOE 

02 

4 

0. 

32000OOOR 

02 

4 

0. 

3  JOOjOOOE 

02 

2 

0. 

34000000E 

02 

4 

0. 

350000 J JE 

02 

3 

0, 

37000000E 

02 

2 

0. 

3MOOOOUOE 

02 

2 

0. 

39O000OOE 

02 

1 

0. 

4  0 JOOOOOF. 

02 

2 

0. 

4  1C0OOOU5 

02 

2 

0. 

430000O0E 

02 

1 

0. 

44000000? 

02 

1 

0. 

45000000'! 

02 

4 

0. 

4o  0 000002 

02 

3 

0. 

47000000E 

02 

1 

0. 

UaOOOOGOE 

02 

2 

0. 

49C000O02 

02 

1 

0. 

5COOOOUOF 

02 

1 

0. 

5  1300OO0E 

02 

2 

0. 

5200COOOE 

02 

1 

0. 

530000002 

02 

1 

c. 

54000000E 

02 

3 

0. 

5500000oE 

02 

2 

0. 

So J00OOOE 

02 

2 

0. 

5H00000OE 

02 

2 

0 

590000  JOE 

02 

1 

0 

62000000E 

02 

■> 

c 

63000OOOE 

02 

3 

0, 

6  4  000  0  00E 

02 

1 

0 

b50000o0E 

02 

2 

0. 

61)0  000  00  E 

02 

1 

c 

6600O00OE 

02 

3 

0 

,e>90G000QE 

02 

1 

0 

7C0OOOO0E 

02 

1 

0 

7300OOOOE 

02 

1 

0 

74000000E 

02 

2 

0 

75000000E 

02 

1 

0 

79000000E 

02 

3 

0. 

83000000E 

02 

1 

0 

84G00000E 

02 

1 

0. 

80O0OOOOE 

02 

1 

0 

. saouooooE 

02 

2 

0 

39000000E 

02 

1 

c 

.  90000000E 

32 

2 

c 

9  1000000S 

02 

1 

c 

. 930OO0O0E 

02 

1 

0 

. 95000000E 

02 

1 

PERCENT 

0.  19048 
0.08036 

PROB&Bn 

0.04167 

0.03274 

»•»»»•» 

O.0253O 

•  *  «** 

0.0  1637 

*•• 

0.  01488 

#•# 

0.01786 

•  •«« 

0.02083 

***• 

0.01339 

*** 

0.01488 

*•* 

0.01637 

•  •» 

0.00893 

** 

0.00893 

•  • 

0.00893 

*  « 

0.  01  190 

*•  • 

0.01  190 

*  ** 

0.  00744 

«  * 

0.01786 

•  »»• 

0.00149 

0.00744 

** 

0.00744 

»• 

0.00446 

* 

0.01042 

•  * 

0.  00446 

* 

0.00446 

• 

0.00298 

* 

0. 30446 

* 

0. 00744 

«• 

0. 00A9  3 

«» 

0.00595 

* 

0.00595 

* 

0.00298 

* 

0.00595 

« 

0. 004  46 

* 

0.00299 

* 

0.00298 

* 

0.00  149 

0. 0079H 

• 

0.0029  8 

« 

0.00149 

0.00149 

0.00595 

* 

0. 00446 

* 

0.00149 

0.00298 

• 

0.00149 

0.00149 

0.00298 

* 

0.00  149 

0.00149 

0.  00446 

* 

0.00298 

* 

0.00299 

» 

0.00298 

* 

0.00  149 

0.0029O 

« 

0.00446 

* 

0.00149 

0.00298 

• 

0.00  149 

0.00446 

» 

0.00149 

0.00149 

0.00149 

0.00298 

• 

0.00149 

0.00446 

• 

0.00  149 

0.00  149 

0.00149 

0.00298 

* 

0.001  49 

0.00298 

* 

0.00149 

0.00149 

0.00149 

. 

Figure     19a    -  USING    LIST    WITH    TELEPHONE    DATA    SET     1 


80 


507 
508 
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511 
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515 
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533 
534 
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536 
537 
538 
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540 
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544 
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546 
547 
548 
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555 
558 
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575 
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580 
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582 
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0. 
0. 

0. 
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0. 

0. 
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0. 
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0. 99000000E 
0.  10600000E 
0.  10900000E 
0.  1  1  100000E 
0.  11200000E 
0.  1  1300000E 
0.  1  1600000E 
0.  1  1700000E 
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0.00149 
0.00149 
0.00149 
00149 
00149 
00149 
00149 
,00149 
,00149 
00149 
00149 
00149 
,00149 
00149 
001  49 
,00298 
0.00149 
0.00149 
0.00446 
0.00149 
,00149 
00149 
00298 
00595 
001  49 
,00149 
,00149 
,00149 
0.00149 
0.00149 
00149 
00149 
00149 
00149 
00149 
00149 
00149 
,00149 
,00149 
00149 
.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.00149 
0.001 49 
0.00298 
00149 
00149 
00149 
00149 
001  49 
00  149 
0.00149 
0.001 49 
0.00  149 
0.00149 


Figure  19b  -    USING  LIST  WITH  TELEPHONE  DATA  SET   1 

(cont. ) 


81 


606 
607 
608 
609 
610 
611 
612 
613 
614 
615 
6  16 
617 
618 
619 
620 
621 
622 
623 
624 
625 
626 
627 
628 
629 
630 
631 
632 
633 
634 
635 
636 
637 
638 
639 
640 
641 
642 
643 
644 
645 
646 
647 
648 
649 
650 
651 
652 
653 
654 
655 
656 
657 
658 
659 
660 
661 
662 
663 
664 
665 
666 
667 
668 
669 
670 
671 
672 


0.  12980000E 
0.  13C50000E 
0. 1328000 0E 
0. 13480000E 
0. 1355000CE 
0. 14120000E 
0.  14290000E 
0.  14890Q00E 
0.  149300005 
0.  15100000E 
0.  15190000E 
0.  154700002 
163300UOE 
173700005 
2C720000E 
243300005 
290600005 
0. 296200005 
0. 3C2600OGE 
0. 3^010000E 
0. 359300005 
0. 368500005 
0. 395200G0E 
0.  4  15700005 
0. 446900005 
0. b20d00005 
0. 76  1400006 
0. 832200005 
0. 9C15000C5 
0. 96250000E 
0. 980600005 
0. 98180000E 
0.  10154000E 
10398000E 
1045100G5 
109390005 
1  12800005 
1  34470005 
143850005 
15155000E 
0.  15294000E 
0. 155040003 
0. 158470005 
0. 158530005 
0. 16280000E 
0- 16299000S 
0. 16361000E 
0. 164 0800 OE 
0. 16817000E 
0. 171740006 
0. 176670005 
0.  182180005 
0.  186490005 
0.  194610005 
0. 218480005 
0. 234990005 
0. 246920005 
0. 264430005 
0. 3C9740005 
0. 356440005 
0. 380030005 
0. 401310005 
0. 47120000E 
0. 475920005 
0. 6  17100006 
0.69775000e 
0- 859930005 


04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
04 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
05 
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V.    ASSESSMENT  CF  VARIABILITY  5CUTINES 


A.    INTRODUCTION 


In  the  previous  sections  Probability  Plotting 
Subroutines  were  presented  for  use  in  informal  estimation  of 
the  form  of  the  distribution  of  a  set  cf  independent 
observations.  In  the  present  section  subroutines  computing 
the  variability  of  some  basic  statistics  are  described. 

We  knew  that  X  (sample  mean)  is  an  unbiased  estimate!  of 

the  population   mean   (    =    x  f  (x)  dx  ) ,   where   X   is   a 

X 

population   with   unknown   distribution  (F  (x)  )  and  unknown 

X 

mean,   and   X  ,  X  ,  ...,  X     is   a   sample   cf   independent 
1    2        n 

observations    from   X.    Furthermore   we   can   compute   the 

Variance  cf  I  as: 

X. 

Var[I]  =  VarC.r-"-] 
1=1  n 

1       n 

= Var[  X     X  ] 

2      i=  1   i 

n 

cl 

=  — -  Var[X]  =  —   . 
2  n 

n 
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Thus   we  see  it  is  possible  to  estimate  the  Var[X],  even 

if   the   population   Variance   (  C  )   is  unknown,  by   simply 

2  2 

using   the   sample   variance   S   as  follows:   Var[Y]  =  S  /n, 

2   rt      _  2  Jsl 

where   S  =  %   (X   -X)   /  (n- 1 )  (an  unbiased  estimator   of  (T~)  • 

2 

Therefore   S  /n   is   an   unbiased   estimator  of  Var[I]  also. 

Furthermore,  especially  for  a  Normal  population,  we  can 
cttain  a  confidence  interval  for  the  nean  using  the 
t-statistic . 


Eut,  although  for  I  there  exists  a  direct  assessment  of 
variability,  for  estimates  of  other  population  parameters 
such  as  Skewness,  Kurtcsis,  Coefficient  of  Variation,  and  so 
forth,  this  is  not  so  simply  obtained.  Thus  several  methods 
have  been  introduced  to  obtain  assessments  of  variatility 
and  for  two  of  them  (Sectioning  of  data  and  the  Jacknife) 
the  Subroutines  'SECTN'  and  'JACK1  will  be  presented. 


E.   SUBfiCuTINE   SECTN 


1 .   Eescri£tion 

'SECTN'  Subroutine  is  used  for  assessing  variatility 
of  estimates  from  data  based  on  the  Sectioning  Method. 

The  tasic  idea  of  this  method  is: 

Assume  we  have  n  independent  observations    X  ,   X  , 

1     2 


•  / 


X     from    a   population   with   unknown   distritution 
n  f\ 

function   F  (x).    Let  a     be   a   parameter   of    F  (x)    and 

8<n)  =  9(X  ,X  ,...,X  ]   be  a  statistic  which  estimates  9  . 


1   2       n 
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Now  we  want  to  estimate  the  variance  of   9(n)    which   is   a 
new  Random  Variable  and  we  work  as  follows: 

1.  Divide  the  sancple  into  r  disjoint  sections,  of 
size  m  (r  should  be  such  that  n=mr.  If  this  is  not  possible 
seme  of  the  latter  data  is  discarded).  ^ 

2.  For    each   section   form   the   same   estimate  C7i  (m)  , 

3.  Compute   the  average  of  9{(.**)#that  is: 

1   r  X 

9  (n)  =  — -  >I  9.  (m)   , 
r   1=1  i 

doing  so,  we  have  an  estimate  of  9  and 

-  1    r  X 

E[§(n)  ]  =  E[— -  SO  (m)  ] 
r   i=1  l 

1     ~ 
=  —  rE[S(m)  ] 

r 

=  E[8>(ro)  ]   , 

also  we  have: 

1   r  "J 
Var(9(n)  ]  =  Var[ Z  9  (m)  ] 

r  i=1  i 

1       r  o 

= Var[  27  9  (m)  ] 

2     i=1  i 

r 

1  Q 

= Var[  y(m)  ] 

r 

= Z    (d.   (m)-tf(n)) 

r   r-1  i=1   l 


1    2 
r    9(m) 
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Thus  we  see  that  Var[  9  (in)  ]  and  Var[  9  (n)  ]  can  be 
unbiasedly  estimated  from  the  sample  variance  of  $i(m)  ,  the 
main  advantage  and  the  purpose  of  the  sectioning  of  data 
method. 

Eur  the  main  disadvantage  of  this  method  is  that   we 

would   like   r   to  be  as  large   as  possible  in  order  tc  make 

2 
the  variability  of  the  variance  estimate   S    as   small   as 

possible.     This    however    may   be   worse   for   the   tias 

properties   cf   the   estimation   procedure.    Therefore   the 

choice  cf  r  is  a  factor  that  should  be  considered. 

It  should  be  noted  here  that  if  (m)  are 
approximately  normal  variates  (This  can  te  tested  ky  the 
previous  described  Subroutine  'NORMPL')  then  confidence 
intervals  for  the  unknown  parameter  can  te  obtained  based 
en  the  t-Statistic,  in  the  following  way: 

5         9(n) 

y(n)  ± t 

1/2      (1-*/2),  (r-1) 

r 

where : 

9  (n)   is   the   mean  cf  the  sectioned  data  statistics 

obtained  from  the  column  named  'mean'  cf   the   second   table 

(•estimated   parameters   of   the   sample  parameters')  cf  the 

program  output. 

1/2 
5  -     /r      is    the   standard   deviation   of   the 
9(n) 

sectioned  data  statistics  divided  by  the  square  root  of   the 

number   of   sections,   obtained  from  the  last  column  of   the 

same  table  of  the  program  output. 

t  is   the   (1-  <*/2)   guantile   cf   the 

(1-  cc/2)   ,  (r-1) 

t-distributicn  with  r-1  degrees  cf  freedom. 
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2.   Program Structure 

'SECTN',  using  the  Sectioning  method,  estimates  the 
following  statistics:  Mean,  Median,  Variance,  Standard 
Eeviaticn,  Coefficient  of  Variation,  Skewness,  Kurtosis, 
Minimum,  and  Maximum.  (The  formulas  which  have  been  used  to 
compute  these  estiaates  are  described  in  the  comments  of  the 
Subroutine.)  The  first  table  is  then  printed  by  the 
program,  containing  the  values  of  these  parameters.  Then 
•SECTN',  using  the  computed  estimates  of  all  sectior.s  for 
each  parameter,  estimates  the  Mean,  Median,  Variance, 
Skewness,  Kurtosis,  and  Standard  Deviation  divided  fcj  the 
square  root  of  the  number  of  sections.  The  second  table  is 
printed  containing  these  values. 

There  are  three  restrictions  in  using  'SECTN': 

1.  The  number  of  sections  must  be  no  greater  than  100. 
If  it  is,  a  diagnostic  messace  is  printed  and  only  estimates 
from  unsecticned  data  will  be  given. 

2.  The  number  of  data  values  must  be  greater  than  3; 
otherwise  a  diagnostic  message  is  printed  without  any 
calculation. 

3.  The  size  of  each  section  must  be  greater  than  2.  If 
it  is  less  than,  or  equal  to  3,  then  the  program  gives 
estimates  for  the  entire  set  of  (unsectioned)  data.  A 
diagnostic  is  printed. 

A  complete  description  of  hew  'SECTN'  operates  is 
given  in  the  Subroutine.  Furthermore,  a  Summary  is  given  by 
typing  en  the  terminal  the  command  DESCRIBE  SECTN  under  the 
CMS  environment.  The  following  response  will  be  given  on 
the  terminal  when  the  user  types  the  above  command: 
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SOBHCUTINE   SECTN 

•SECTN*  is  intended  to  estimate  a  set  of  tasic 
statistics  cf  a  given  set  cf  observations  using  the 
•Sectioning  cf  data  method'.  Also  for  each  estimated 
statistic,  estimates  of  some  basic  statistics  such  as  the 
mean,  standard  Deviation,  and  so  forth,  are  given. 

It  is  called  by  : 
CALL    SECTN   (  X,  N,  K  ) 

where : 

X    Is  the  array  of  data 

N    Is   the   nunsber   cf   data  values  (must  be  greater 

than  3) 
K    Is  the  number  of  desired  sections  (no  greater 

than  100) 

Note:  k  should  be  a  number  which  minimizes  the  number  of 
data  points  that  will  have  to  be  discarded.  'SECTN'  places 
the  data  into  the  equal  size  sections  discarding  any  data 
left  over. 

Fcr  k<3  or  k>1C0  or  (n/k)<3  only  estimates  from 
unsectioned  data  will  be  given  and  no  estimates  fcr  the 
estimated  statistics  will  be  computed. 

No  output  is  expected  if  n<3  . 

Kcre  information  is  given  in  the  subroutine. 
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3.   Using  1SECTNJ.  with  Telephone  Data  Set  J 

•SECIiN'  was  used  on  Telephone  Data  Set  1  to  assess 
the  variability  in  the  Mean,  Median,  Variance,  Standard 
Deviation,  Coefficient  of  Variation,  Skewness  and   Kurtcsis. 

The  672  data-points  cf  Telephone  Data  Set  1  were 
broken  dcwn  into  16  disjoint  sections  with  42  data-points 
per  section.  Because  of  this  break-down  no  data-points  were 
discarded.  Compare  the  values  cf  row  ' unsectioned '  (see 
figure  21)  with  the  values  of  the  corresponding  statistics 
computed  by  the  Subroutine  ' HISTGS/HISIFS •  (see  figure  1)  . 
The  values  are  the  same. 

Now  if  we  want  to  assess  the  variability  in  any  of 
the  aoove  parameters  we  proceed  as  follows: 

1.  He  take  the  mean  (0  )  of  the  parameter,  whose  the 
variability  we  want  to  a.ssess,  from  the  table  under 
•ESTIMATED   EAEAMETEfiS   OF  SAMPLE   PAHAMETEES'. 

2.  We  rake,  from  the  same  table,  the  value  (S)  frcm  the 
last  column  (STD.  DEV./NS**.5)  for  the  same  parameter. 

3.  Using   the  t-statistic  with  k- 1  d.f.   and  the  formula: 

B  ±  St 

(1-^2)  ,  (k-1) 

we  get  a  (1-0-)%   confidence  interval  for  the  parameter. 

Example 

In   order   to  find  a  95X  confidence  interval  for  the 

Skewness  of  Telephone  Eata  1  we  take: 

8=  4.9343,     S  =  .30312,     t         =  2.131 

.975,  15 

Therefore    the    95%    confidence    interval    is: 

4.9343    ±    .30312*2.131    =    4.9343    ±     .64595 
===>    [ 4.2884,     5.58C2]. 
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It  should  be  ncted  here  that  the  use  of  the  Variance 
estimate  from  the  sectioned  data  in  order  tc  get  confidence 
intervals  of  the  parameters  is  based  en  the  normality 
assumption  and  the  independence  of  the  estimates  frcm  the 
sections.  The  normality  will  depend  on  the  number  of 
data-points  in  each  section,  which  should  be  kept  large. 
This  reguirenent,  however,  conflicts  with  the  need  to  make 
the  number  of  sections  large  to  reduce  the  variability  in 
the  estimate  of  the  variance  of  the  statistics.  The 
skewness  estimates  of  each  of  the  16  sections,  can  fce  run 
through  the  normal  plotting  routine  to  see  whether  the  use 
of  the  t-statistic  confidence  interval  is  valid. 
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C.   SUBROUTINE   JACK 


1  •   Description 

This  routine  is  used  for  assessing  variability  of 
estimates  of  data  based  on  the  Jacknife  method,  and  for 
reducing  the  bias  in  estimates  also.  It  is  particularly 
useful  fcr  data  with  small  sample  size. 

A  big  picture  cf  the  method  has  as  fellows: 

Let  X  ,  X  ,  ...,  X   be  a  sample   of   n   independent. 
1    2        n 

identically   distributed  observations  from  a  population  with 

unknown  distribution   function   F  (x) .    Alsc   let  u    he   an 

X 

unknown   parameter   of  F  (x)   tc   be  estimated.  Furthermore 

X 

suppose  a  method   (biased   or   unbiased)   is   available   for 

estimating  u   .  Then  we  proceed  as  follows: 

1.  Divide  the  sample  size  into  r  disjoint  groups,  each  of 
size  m  (r  should  be  such  that  n  =  ir.  If  this  is  not 
possible  seme  of  the  last  data-points  will  he  discarded.) 

2.  Compute  the  estimator  8  (n)  of  v)  based  on  all  n  =  mr 
observations.  ^ 

3.  Compute    the   estimator  B  (n-m)   based   on   the   n-m 

i 

th 
observations ,  having  deleted  the  i   group. 

A. 

4.  Compute  the  so  called  'PSEUDO   VALUES'  d    (n) : 

i 

9  (n)  =  r  9  (n)  -  (r-1)  9#  (n-m)  ,     i=  1,2,. ..,r  . 
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Then   the  Jacknife  estimator  is  defined  to  fce  the  average  of 
Eseudo  Values: 

a       1    r  A 

9<n)  =  —   2  8  (n) 
r    i=1  i 


1     2^ 
— -  (  r  0(a) 


(r-1)  2f  9  (r-m)  ) 
i=1  i 


q       r-1 
r  c7(n) 

r 


X     B    (n-m) 
i=1   i 


Havicg   so   defined  the  Jacknife  estimator  it  can  be 

proved  that   this  estimator  is  an  unbiased  estimator   cf  &, 

-2 
except   for   terms   n    and  higher  order,  assuming  that  the 

bias  has  the  form  cf  E[  9(n)  ]  =  0  +  a/n  +   terms   of   higher 

crder   in   n.    That  is,  the  Jacknife  estimatcr  eliminates  a 

n    bias   term.  Namely  if 

E[  9  (n)  ]  =  O  +  an    +  0  (n   )   , 

then  we  have: 

A  ^        r_-j    r  ~ 

E[9(n)  ]  =  E[r$(n) X    $  (n-m)  ] 

r    i=1  i 


=  rE[8(n)  J  - 


C"1    r   ^ 

X  E[9  (n-m)  ] 

r    i=1    i 


=  r  ( 


-2 


r-1 

r 


-2 


q     -1      -2      r-l   q.   ra     r     -  i. 
C7  +  an   +  0  <n   )  ) (rtf+ + C  (n   )) 


n-a    n-m 


^    ra  q      a 

ra+ +  ...-(r-1)  (  0  + +...) 

n  m  (r-1) 


=  r 


n    ra         q    (r-1)  a    Q 
n  m(r-1) 
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q    q    ra     a    q 

lU  -   ru   + +  O   ♦ 

mr     is 


=  9. 


A  special  case  of  Jacknife  estinator  can  be 
mentioned  here,  called  the  'Complete  Jacknife  Estimator', 
where  r=n  and  m-1. 

Seme  properties  of  the  Jacknife  estimator  foilcw: 

-2 

1.  It  is  ucbiased  up  to  order  n 

2.  The   ' Eseudo  values'   can   be   used  to  obtain  variance 

estimates  for  the   Jacknife   estimator   since   they   can   be 

considered   (see   [18])    as   approximately   independent  and 

-1  r  a     *•         2 
identically   distributed.   Thus   (r(r-1))     ST  (9  (n)-9(n)) 

i=1   i 

A. 

should  be  an  approximate  estimate  of  Vai[9  (n)  J  , 

and 

*  1/2 

i    l  S      2  1/2 

£  2.    {&.  (n)  -  V(n))    } 

i=1   l 

should  be  approximat elly  distributed  as  a  t-statistic  with 
r-1  d.f.  This  procedure  is  particularly  useful  if  the 
number  n  cf  data-points  is  small,  but  it  must  be  used  with 
care. 

3.  For  large  n  it  can  be   shown   that   for   the   ccnplete 
Jacknife  estimator,  under  very  general  conditions  we  have: 

Var[  8<n)]  >   Var[  D  (n)  ]  . 
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2.   112915™ Structure 

Subroutine   'JACK'   is  a  FORIRAN-callable  Subroutine 

which  takes  a  set  of  data,  groups  it  into  r  disjoint   groups 

th 
and  for  each  set  of  the  r-i    group  computes  and  prints  in  a 

table  the   statistics;   Mean,   Median,   Variance,   Standard 

Deviation,   Coefficient  of  Variation,  Skewness  and  Kurtcsis. 

The  same  statistics  are  also  computed  for  the  ungrouped   set 

of    data.    Then,   using   these   statistics,   the   Jacknife 

estimator  of  the  above  parameters  are  computed   and   printed 

in   anotner   table   along   with   its   Variance   and  Standard 

Deviation . 

The  program  is  divided  into  the  main  program  ar.d  tne 
Subroutine  JACKES.  The  main  program  groups  the  data  and 
successively  calls  the  Subroutine  'JACKES'  in  order  to  get 
the  estimates  of  the  above  parameters  for  ungrouped  data  and 
for  each  group  as  well.  A  table  is  then  printed  containing 
the  estimated  parameters.  The  'Pseudo-Values'  for  each 
parameter  are  also  computed  by  the  main  program,  and  are 
used  to  compute  the  Jacknife  estimator  for  each  parameter 
and  its  Variance  and  Standard  Deviation. 

To  avoid  division  by  zero  the  number  of  observations 
must  be  greater  than  three,  otherwise  no  output  is  given. 
For  the  same  reason  the  expression  (r-1)  (n/r)  must  evaluate 
to  greater  than  three,  otherwise  the  program  will  give  only 
estimates  of  ungrouped  data. 

A  complete  description  of  how  'JACK'  operates  is 
given  in  the  Subroutine.  A  Summary  is  printed  on  the 
terminal  by  typing  DESCRIBE  JACK  under  the  CMS  environment, 
which  responds: 
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SUBROUTINE   JACK 

Jack  is  intended  to  estimate  the  statistics  :  Mean, 
Median,  Variance,  Standard  Deviation,  Coefficient  of 
Variation,  SJcewness  and  Kurtosis  of  a  given  set  of 
independent  cbser vaticns,  using  the  Jacknife  method.  In 
addition  to  the  above  parameters  which  are  given  fcr  each 
group,  the  Jacknife  estimator  along  with  its  variance  and 
standard  deviation  is  given  fcr  each  parameter. 

It  is  called  by: 

CALL    JACK   (  X,  XS,  STAT,  N,  IG  ) 
where: 

X     is  the  array  of  data  of  dimension  N 

XS    is  a  work  array  of  dimension  N  (returns  ordered 

data) 
STAT   is  a  two-dimensional  work  array,  dimensioned  by 

<IG,7) 
N      is  the  number  of  data-values 
IG    is  the  number  of  groups  plus  one  (=r+1) 

Note  :  r  must  be  such  a  number  as  to  minimize  the 
number  of  data-points  that  will  have  to  be  discarded.  JACK 
places  the  data  into  the  equal  size  groups  discarding  any 
data  left  ever  (last  observations  in  order  of  original 
input) . 

Nc  output  is  expected  if  N  <  3.  Furthermore  if 
IG  <  2  or  (IG-2)  (N/(IG-1))  <  3  only  estimates  for  ungrouped 
data  will  be  given  . 

Wore  information  is  given  in  the  Subroutine. 
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3.   Qsing  iJACKJ.  with  Telephone  Data  1 

Tc  assess  the  variability  in  the  Mean ,  Median, 
Variance,  Standard  Deviation,  Coefficient  of  Variation, 
Skewness  and  Kurtosis  cf  Telepfone  Data  1  the  Subroutine 
•JACK'  is  now  used.  In  order  tc  avoid  to  discarding  any  of 
the  data-points  we  use  16  groups  with  42  data-points  per 
group . 

Observing  figure  22  we  see  that  the  values  of  the 
parameters  for  ungrouped  data  are  exactly  the  same  as  the 
corresponding  values  of  the  same  parameters  of  figure  1 
which  was  produced  by  using  the  Subroutine  HISTGS/HISTFS. 
Now  in  order  to  assess  the  variability  in  anv.  of  the  above 
parameters  we  have  to  use  the  t-statistic  with  (r-1) 
d.f.  along  with  the  values  printed  under 
•ESTIMATED  JACKNIFE  PARAMETERS'  using: 

9  ±  st 

(1-^2)  ,  (r-1) 

where  C7  is   the   Jacknife   estimator   and  S  is  its  Standard 

Deviation.  Fcr  example  to   assess   the   variability   of   the 

Skewness  with  a  confidence  level  €K=    .05  we  have: 

&  =  7.37321,      S  =  .900118,     t         =  2.131  , 

.975,15 

Therefore  a  S5X  confidence  interval  of  the  Skewness  is: 

7.37321  ±  .900118*2.131  =  7.37321  ±  1.S1815 

===>   [5.4551,  9.2914]   . 

This  should  be  compared  with  the  point  estimate  4.9343  and 
the  confidence  interval  estimate  [4.2884,  5.5802]  obtained 
from  subroutine  SECTN.  The  data  is  so  skewed  that  one  will 
inevitably  have  trouble  here  with  the  Jacknife  procedure. 
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Figure    22    -         USING    JACK    WITH    TELEPHONE    DATA    SET    1. 
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VI.    SINGLE^ EBVEE  FISS INCOME  FI3ST-SZHVED  QUEUES  fllTH 

•EARM&«  STRUCTURE 


A.    THE  EAEKA<P,Q)  PROCESS 


The  starting  point  for  these  processes  is  the  definition 

cf  a   first-crder   autcregressive   model   for   a   stationary 

seguence  cf  random  variables  {X  } : 

i 


X  =  rX    +  8  ; 
i     i-1     i 


i  =  0,±1 ,±2,  . . 


If  the  marginal  distritution  of   the   X     is   fixed   to   be 

i 
exponential  with  parameter  7\   for  all  i; 


->x 
P  {X   <  x}  =  1-e    , 

i 


>>  0,  x>  0   , 


then  8   is  zero  with  probability  r  and   exponential {y)       with 
i 

prooabiiity  T**r.   Thus 


rrX 


i-1 


^rX     +  E 
i-1    i 


w.  p.   r 
w.  p.   1-r , 


(D 


where  {E  }  is  a  seguence  of  i.i.d.  exponential  (DO  random 
variables.  The  process  defined  by  (1)  is  the  exponential 
autcregressive  process  of  order  1,  the  EAR(1)  process.  The 
correlation  structure  is  r(j)  =  r  ,    j  =  0,1,2,...  . 
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The  first   order   moving   average   expcrectial   process, 
EMA(1)  is  defined  as 


i 


4e  v. p.  & 

i 

il       +  E       w.p.  (1-#) 
i    i-1 


(0 


<£   <  1;   i  =  0,±1,±2,...) 


and  has  correlations  o(1)  =   (1-$)  and  p(1)  =  0,  j  =  2,3,... 


The  general  E2IA  (g)  model  takes  the  form 


q  i 

g  i    ug-1  i-1 


X   =  S 


Cq  1    ^q-1  i-1 


+  ^E 

1  i-q+1 


w.p.   b 


/t   ♦  /   E     +...♦©  E       +  E       w.p 
g  i    g-1  i-1  1  i-g+1     i-g 


for  0  <  ^>    ,  £  ,  ...,  tf    <    1;   i  =0,12,  ...   wher« 
1    2         q 


b   =  S 

i 


4 

(1-^)  ...  (1-4)4 
g  ii-1 

(1-^)  ...  (1-^) 

q       i 


g+1 


w.p.   b 


(2) 


w.p.   b 


i  -  g+1 , 
g  >  i  >  2, 
i  =  1   . 


The  serial  correlations  are  given  by  the  equation 


P 


(q) 


(j)  =  corr  (X.  ,X.   )  =  \ 

l   l-r 


q-j+i  ,  u 

X  b  b 

V=j  V  v+j 

0 


(1  <  j  <  q) 


(g+1  <  j  <oo) 
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Tc   convert    this   EMA  (g)    process    into    a    mixed 

autoregr essive   moving-average  exponential  process  of  orders 

1  and  g,  called  EMA(1,g)  we  replace  E    in  (2)  with 

i-g 


rA 


A     =  • 
i-g 


i-g-1 


rA       +  E 
i-q-1     i-q 


w.p 


w.p.   (1-r) 


For    further   results   en   these   processes   and  their 

properties  see  Gaver  and  Lewis  (1978),   Lawrance   and  lewis 

(1977),   Jacobs   and   Lewis   (1977),   and  Lawrance  and  Lewis 
(1978)  . 


£.    USE  CE  AUTOREGBESSIVE  (EAR(P)),  MOVING  AVERAGE  (EKA(Q)) 
AND  MIXEE  ST5UCT0BES  (EAfiMA(P,QJ)  IN  MODELLING  QUEUES 


Consider   for   simplicity   a   queue   with  a  single  input 

stream  and  a  single  server,   and   a   f irst-ccme-f irst-served 

(EIEO)   service   discipline.   Let  S  ,  i  =  0,1,2,...  ,  denote 

i 
th 
the   service   time   for   the   i    arrival,   and    let    X  , 

i 
th 
i  =  1,2,...  ,   denote   times   between  arrival  of  the  i    and 

st 
(i-1)    customers.   As  is  usual  we   assume   that   the   first 

customer   (with   service   time   S  )  arrives  at  time  zero  and 

0 

finds    the   gueue   empty. 

If      the       {S   }       and    {X    }     sequences    are    i.i.d.    exponential 
i  i 

random    variatles      with      parameters   ^  and  a   respectively,      we 
have    the    H/M/1    queue. 
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Now  let 


E    be  exponential (>)  and  independent,  i  =  0,+1,...  ; 
i 

£   be  exponential (a)  and  independent,  i  =  0,+1,...  . 
i 


We  want  to  model  queues  with  correlated  (autocorrelated 
and/or  cross-correlated)  service  and  inter-arrival  times, 
the  service  and  inter-arrival  times  both  having  marginally 
exponential  distributions.  We  also  want  the  gueing  model  to 
include  the  a/M/1  gueue  as  a  special  case. 

There  are  five  simple  possibilities  which  we  pur  forward 
here,  based  on  the  use  of  EARMA(p,g)  processes,  giving  what 
we  call  an  EAfiMA  MD/MD/I  queue. 


In   what   follows   let   p    and  p   be,  respectively,  the 

S       X 

order  of  the  autoreqressive  components  of  the   service   time 

sequence   {S  }  and  the  inter-arrival  time  sequence  {X  }  ,  and 
i  i 

let  q  ,g   be,  respectively,  the  order  of  the  moving   average 
S   X 

component   of   {S  }   and   {X  } .    These   parameters  can  take 

i  i 

values   C,1,...  .    If   p  =  0,   g  =  0,    the    sequence    is 

independent;   if  p  =  0 ,  the  process  is  purely  moving  average 
and  if  g  =  0,  then  the  process  is  purely  autcregresive . 


The  five  possibilities  or  cases  are  as  fellows: 

1)   Let  {S    }  be  EARMA(p  ,g  )  over  {E  ,  2    ,  ...}; 
i  S   S  i    i-1 

let  iX.}  =  {£.}. 
l      l 

Thus  tne  arrivals  are  a  Poisson  process  and  the  service 
times  are  autoccrr elated,  i.e.  a  mixed  autoregr essive 
moving  average  exponential  sequence. 


106 


2)   Let  S  ■  E    ,       i  =  0,1  ,2,  ..  .  ; 
i    i 

let  X  fc€  EARHA(p  rq  )  ever  {£  ,  £   ,  ...} 
i  XX  i    i-1 


Then  the  service  times  are  independent  and  the  arrival 
process  is  non-Poisson  because  of  the  dependency  between 
the  inter-arrival  times,  i.e.  the  arrival  process  is  a 
point  process  with  EARMA(p,q)  structures. 


3)   Let  (S  }  he  EARMA(p  ,q  )  ever  {E  ,  E    ,  ...}. 
i  S   s         i   i-1 

Let  {X  }  te  EARMA(p  ,q  )  ever  {  %   ,  g        ,  ...}. 
i  XX         i   i-1 

The  service  and  arrival   processes   are   autocorrelated, 

but   the   two   processes   are  independent.   The  marginal 

distributions  of  {S  }  and  {X  }  are  still  exponential. 

i        i 


4)  To  ccuple  the  two  processes,  with  resultant  dependence 
in  the  arrival  and  service  processes,  the  sinplest 
procedure  seems  to  be  the  following. 


Let  {S  }  be  EARMA(p  ,q  )  in  the  following  sense: 
i  S   S 

{S  }  is  EASMA(p  ,q  )  over  {E  ,  S.&,  a.  fc   ,  ...}  • 
i  S   S         i   y    i   ^  i-1 


Then  if  X  =  t  ,      i  =  0,1,2, we  have  that  (S  }  is  an 

i     i  i 

autcccrr elated  sequence  and  also   cross-ccrrelated   with 

{X   =  6  },    i.e.    (S  ,X  }   is   a   bivariate   dependent 
i    i  i   i 

sequence  cf  random  variables  with   exponential   marginal 

distributions,   although   {X  }   is   still  Poisson.   Note 

i 

that  the   S    sequence   is   an   autccorrelated   sequence 
i 

because    of    the   cross-coupling,    but   not   a   pure 
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Let 


EARMA(p    ,q    )    sequence.         It      does      have      an      exponential 
S      S 

marginal   distribution    however. 


As    an   example    let    p      =    1,    q      =    1. 

S  S 


S      =    E 

0  0 

S      =    ft E                       w.p.      /  A      =  a  & 

1  S     1                                           S  1^1 

=    d  E  +    A  w.p.        (1-^  ) 

S^    =    ft  E                       w.p.       /f            A  =    rA  w.p.    r 

2  S    2                                      ttS              2  1 

=   ^  E  +    A            w.p.        (1-/)  =    rA      +  £L£        w.p.      (1-r) 

S    2           2                                   S  1         ^    2 


=    #    E  w.p./ 
Si  s 

=    /f  E      +    A  w.p.       (1-      )  =    rA  +a&      w.p.      (1-r) 

^Si  i  S  i-1       ^    i 


S      =    16   E  w.p.      £  A      =   rA  w.p.    r 

i  S    i  S  i  i-1 


In   these    equations    one    can    replace    6      with   X        and      then 

i  i 

the      cross-coupling      between   the    sequences    {S    }    and    {X   } 

i        i 

becomes  apparent. 


Interpretation.   We  have  positive  correlation  between  S 

i 

and      q        previous    inter-arrival   times     (when    p      =    0) .      If 
S  S 
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the  6  s  (j  =  1,1-1,.  . . ,    i-p  -1) 


are   short,   then   S 


will  fce  short  if  all  the  g   previous  inter-arrival  times 

S 

are  short.  This  models  the  case  where  the   server   tends 

to   i_peed   U£  if  the  c[ueue  gets  lo_ng.   Cf  cource  h€  also 

slows  dcwn  when  the  queue   gets   shcrt   and   it   is   not 

immediately   clear  what  the  effect  on  an  average  waiting 

time  will  he. 


4«)     Similar         to         4         but  {X    }         is         EABMA(p    ,q    )         over 

i  XX 

£  ,  2.2  ,   2LE    '•••   and  S   =E. 
i   c*  i   ex  i-1  i    i 

Interpretation  is  net  clear,  but  it  would  have   the 

same   affect  as  balking  in  the  input  stream.   If  service 

times  get  long  (and  presumably  the  queue  gets  long)  then 

inter-arrival  times  get  leng. 

5)  A  mere  general  possibility  which  allcws  one  to  model 
dependency  in  the  input  stream,  in  the  service  stream 
and  dependency  which  couples  them  is  as  fellows: 


Let   {E  }   be  a  unit  exponential  independent  sequence,  a 
i 

driving  sequence. 


Let     {1   }  be    EARMA(p    ,q    }     over     {E    ,    E         ,...}. 
i  S      S  i        i- 1 

Let     {1   }  be   EARMA(p   ,g   }    ever    {£  ,     p.       #•••}• 
i  XX  i        i-1 

Let     {3    }  be    EARMA<p      ,q      )    over    {2    ,    ^D    ,    ^D         ,...} 
i  ci       cL  i  i  l- 1 

Let     {X    }  te    EARMA<p      ,q      )     over    {I    ,    W,,    aD,        /...} 
i  ct      c£  i  l  i-1 
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This   is   a   very  general  scheme  which  reduces  tc  all  of 
the  previous  schemes  as  special  cases.   For   example   if 

p    =  g    =  p    =  g    =C,   then   we   have   case   3.   If 

Cl       CL       Gi       ci 

g    =1/  f    =  0/   /   =  Or  P   =  g   =  0;  p   =  0  (and   p 
cL  cz  ct  S     S        X  ct 

and  g    are  not  egual  to  zero)  we  have  case  4. 
CL 

In  all  cases  the  basic  eguation  for  M/M/1  gueue   waiting 

time  W   still  holds: 
n 


W   =  0 
0 

+ 
*  =  ( »  +  £  -  X   )       n  =  1 , 2 , . . .  . 

n+1      n    n     n  +  1 

t 

These   eguations  can  be  used  to  generate  successive  K   s 

n 

in  a  simulation.   We  have  not  touched  on  hew  the   correlated 

sequences  are  started.    This  is  somewhat  artitrary  as  it  is 

not  known  (except  for  the  EMA(1)  and  MA(1)  processes)  hew  to 

start   them   so   as   tc   produce   a   stationary   exponential 

sequence.   The  problem  is  aggravated,  as  in  cases  4,  4    and 

5,  when  we  are  dealing  with  tivariate  exponential  sequences. 
Then  the  aarginal  processes  may  be  stationary  but  not  the 
tivariate  process. 
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C.   PRCGEAM  STEOCTURE 


1 .   General 

Tc  simulate  the  EARMA(p,q)  model,  a  FORTRAN  Picgram 
has  been  written  tc  take  care  of  all  the  existing  cases  of 
that  model  and  also  of  the  H/M/1  queue.  Eepending  en  the 
input  values  which  the  parameters  (as  described  below)  will 
take,  the  program  can  be  used  for  creating  in  each  run  a 
particular  case  of  the  EARMA(p,q)  queue.  Eecause  of  that 
generality  cf  the  program  it  is  not  suggested  that  it  be 
used  for  these  EARMA  cases  where  either  the  arrival  tines  or 
the  service  times  cr  both  are  not  autccorrelat ed  or 
cross-correlated.  The  reason  for  this  suggestior.  is 
efficiency.  It  is  suggested  for  use  with  the  EARM2(p,q) 
cases  where  both  arrival  and  service  times  are 
autccorrelat€d  and,  furthermore,  cross-correlated  with  any 
order.  The  program  is  divided  into  the  main  progran  and 
into  the  Subroutines  BETAS,  AUTOR  and  EARtfA. 

2«   Ma.in  Program 

The   main   program   performs   all  I/C  operatiors  and 
calculations    required    for    the    desired    statistics. 
Specifically  it  is  designed: 
a)  To  read 

1.  All  random  number  generator  seeds  required  to 
generate  the  sequences  of  the  Exponential  and  Uniform 
variates  to  be  used  during  the  execution. 

2.  The  values  for  the  variables: 

N  :   Number  of  arrivals  to  be  generated. 
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a    :   Sumter  of  replications. 

RX,  RS  :   Arrival  rate,  service  rate. 

C£X,  CFS,  CPG  :   Taking  the  values  1  cr  0   (1  £cr  the 

coupled  process) . 
KX1,  £Xy  KS1,  K£  :   The  order  of  moving   average   for 

arrival  and  service  processes. 
QX1,  QX,  QS1,  Q£  :   Taking  the  values  1  or  0    (1  for 

autor egressive  case) . 
EX,  BX1,  BS,  BS1  :   Arrays  with  the  values  of  -qs>  for 

mcving  average  part. 
RH0X1,  RHGX,  RHCS1 ,  RHCS  :   Parameter  values  for  the 

autor egressive  farts. 

t)  lo  calculate  and  print 

SUMX  (I,J)  ;  SUiSXM  (I, J)  ; 

SUtfS  (I, J)  ;  SUMSM  (I, J)  ; 

W  (I, J)  ;  WM(I,J)  ; 

HE  (I, J)  ;  WME  (I, J)  ; 

D  (I, J)  ;  EM(I,J)  ; 

for  I  =  1,2,3,4,     J  =  1,2, ...,H, 

where : 

(In  what  follows  k  gets  the  values  N/4,  N/2,  3N/4,  N 

for  1=1,2,3,4  respectively  and  J=1,2,...,M) 

k   (J)        (J)       th 
SUMX(I,J)  =  Z  X    ,     X     the  r   arrival   for 
r=1  r         r 

th 
the  J    replication  in  tae  EARMA  gueue; 

k      (J)        (J)       th 
SUK2H(I,J)  =  Z  (XM)    ,   (Xil)     the  r     arrival 
r=1     r         r 

th 
for  the  J    replication  in  H/tt/1  Queue; 

k   (J)        (J)   L    th 
SUMS  (I, J)  =  51  S    ,     S     the  r    service  for 
r  =  1  r         r 

th 
the  J    replication  in  the  EARMA  gueue; 
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JL  (J)  (J)  th 

saMSM(I#J)     =    Z    (SB)         ,       (SM)  the   r  service 

r=1  r  r 

th 
for    the    J        replication   in    M/M/1    Queue; 

(J)  (J)  (J)  (J) 

H  (I, J)     =    W  =    max{    (W  +    S  -    X         )  ,    0} 

k  k-1  k-1  k 

th  th 

=  waiting  time  of  k   arrival  in  J   replication 

of  the  correlated  queue; 

(J) 
WM  (IrJ)  =  (WM)       The    same  as     W(I,J)     but 

k 

fcr    a/M/1    queue; 

k       (J) 

wb  (I, J)    =      rw        /k 
r=1    r 

A  (J) 

HHE(I,J)     =       X     (WM)  /k 

r=1  r 

D  (I,J)  =  {S0MS(I,J)  -  SCUX  (I,J)  }/K  -  1/RS  +  1/SX 

DM|I,J)  =  (S0MSH(I,J)  -  SUMXM  (I, J) }/K  -  1/RS  +  1/RX 

k       (J) 
S0flSM(I,J)  =  X     (SM)     , 
r=1     r 


Ncte:  Between  the  group- variables  (KX1,  KX,  KS1,  KS)  ,  (C.X1, 
£X,  QS1,  QS) ,  (CPX,  CPS,  CFG),  (RHOX1,  flfcCX,  RHOS1,  RHCS) 
and  (BX1,  EX,  BS1,  3S)  there  are  the  following  relations: 

1.  The  number  of  elements  of  each  variable  (array)  from 
E-group  should  be  the  same  as  the  value  cf  the  corr espending 
variable  frca  K-grcup. 

2.  The  '0'  value  cf  a  variable  from  Q-grcup  implies  that 
the  corresponding  variable  from  R-group  is  net  needed. 

3.  The  value  ' 1'  of  CPG  dominates  any  value  of  CPX,  or 
CPS. 

Any  valid  combination  in  the  values  of  the  above 
variables  specially  cf  the  C~group,  C-grcup  and  K-grcup  (0 
cr  not  0),  generates  a  particular  EARMA(p,q)  case  and  is 
accepted  by  the  program. 
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The  tasic  approach  in  the  main  program  is  tc 
generate  at  once  and  for  each  replication  three  independent 
exponential  sequences  with  parameters  XL  (arrival  rate),  SI 
(service  rate)  and  1  (unit  exponential)  which  are  stored  in 
the  arrays  EXPXL,  EXFS1  and  EXP1  respectively.  Ther  the 
Subroutine  AUTOR  which  is  involved  with  the  autoregr essive 
part  of  the  model  is  called  (if  it  is  reeded  for  the 
particular  EABMA  case)  .  A  loop  follows  which  is  executed  as 
many  times  as  the  value  of  N.  From  the  loop  the  EAEMA 
Subroutine  is  called  which  is  involved  with  the  acving 
average  part  Gf  the  model.  The  number  of  calls  of  Subroutine 
EARMA  depends  on  the  particular  case  which  is  being  rur.  The 
statistics  and  all  information  that  will  help  us  to  analyze 
the  model  are  also  computed  inside  the  lccp  and  stored  in 
tne  arrays  described  above.  The  program  continues  execution 
until  the  desired  statistics  for  all  a  replications  have 
been  calculated  and  gathered.  Then  the  output  part  of  the 
program  fellows  which  gives  us  the  values  of  calculated 
statistics  on  paper  and  on  punched  cards  fcr  analysis. 

3.   Subroutine  BETAS 


This  is  a  simple  Subroutine  and  its   purpose   is   tc 

return    an    array   SUH3X1/X/S1/S   with   element-values   as 

fellows: 

I 
SUMBX(I)  =  X     EX(r)  ,     I  =  1 , 2 , . .  .  ,  value  of  K-grcup 


r=1 


to  be  used  for  rhe  choice  of  the  order  in  the  moving 
average.  It  is  called  at  most  4  times  during  the  execution 
and  only  if  the  K-group  variables  have  value  greater  than 
•0*. 
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*•   Subroutine  AOTCR 

'AOTCR'  accepts  in  each  call  an  exponential  sequence 

(generated  previously)  of  variates  and  transforms  it  into  an 

autocorr elated  sequence.  It  is  called  only   if   any   of   the 

Q-group  variables  eguals  1.  For  each  Q-group  variable  having 

value  1  it  is  called  M  times. 

5.   2ub£CUtine  EAR  MA 

This  is  the  main  Subroutine  of  the   program   and   it 

has   been   created  so  that  it  can  be  used  for  all  EARMA(p,q) 

cases.   Because  of  the  generality  guite  a  few  parameters  are 

transferred   and   in   each   call  just  one  value  is  returned. 

th  th 

This  value  is  the  I    arrival  or  the  I    service  time  as   it 

has  been  modified  because  of  the  loving  average.  It  is 
called  only  if  the  particular  EARMA(p,g)  case  reguires 
moving  average;  that  is  only  if  any  of  the  K-group  variables 
is  greater  than  0. 

6  •  1  il^  requirements 

Seme  CPU-times  have  been  gathered  in  the  course  of 
the  simulation  of  some  EARMA(p,g)  cases,  using  the  IBM 
360/67,  FORTRAN  H  compiler  as  follows: 


TIME 

9  Mia.  51  Sec. 

22  Min.  25  Sec. 

43  Bin.  30  Sec. 

55  Kin.  20  Sec. 


M 

M 

K- 

-GROUP 

QS1 

C-GROUP 

2000 

500 

0 

1 

0 

500C 

500 

0 

1 

0 

10000 

500 

0 

1 

0 

1000C 

500 

0 

1 

CPS=1 
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from  these  results  vie  see  that  the  value  cf  N  is  the  main 
factor  that  affects  the  CPU-time  (linearly)  .  We  see  also 
that  one  cf  the  C-group  variables  may  increase  the  CPC-time 
by  2555.  We  ignore  the  factor  of  M  (it  is  another  main 
factor  as  N  is)  because  we  may  step  the  progran  at  any 
replication,  getting  the  last  random  number  generator  seeds 
and  then  continuing  another  time  using  these  seeds. 


D. 


ANALYSIS  OF  SIMULATION  EESDLTS 


Since  nc  analytical  properties  of  the  EABHA(p,q)  models 
can  be  derived,  a  simulation  has  been  done  to  study  their 
properties.  Eut  as  we  have  noticed,  this  program  is  a 
general  program  accomodating  all  EABMA(p,g)  cases.  Since 
oar  study  is  limited  to  only  the  two  cases  telow,  in  crder 
to  save  CPU-time  a  modified  program  has  been  usee  that 
simulates  these  particular  two  cases. 

1 .   0.ueue  with  Dependent  Service 

The  model: 


U.}  =  6. 

i  i 

[S  }  is  EARMA  (1 ,0; 
i 


That  is 


o   — 

i 


rA 


i-1 


t  rA     +  E 
i-1     i 


over    {  E  ,  E    ,  .  .  .  } 
i    i-1 


w.  p.   r 


w.p.   (1-r) 


where   A   =  E 
0     0 


To  simulate  this  model  we  chose  the  values  .25,  .50, 
.95,   .99   for   traffic  intensity  t  and  the  values  .25,  .50, 
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.90,  .95,  .96  for  the  correlation  r,  in  order  to  cover  a 
representative  range.  In  what  follows  the  string  S#tt#rr 
stands  for  the  run  with  t=tt  and  r=rr .  That  is  the  run 
S#25#50  stands  for  the  run  with  t=.25  and  r=.50  For  each  t 
and  r  we  ran  the  program  and  the  required  statistics  were 
gathered,  analyzed  and  plotted  using  the  Subroutines 
HISTGS/FS,  NCRMPL  and  EXPLT.  The  sample  statistics  from 
HISTGS/FS  were  also  tabulated  separately  frca  the  figures  so 
that  one  could  obtain  an  overall  picture  of  whether  or  not 
the  waiting  times  had  converged  in  mean  and  in  distrirution 
to  tne  Uniting  distritution . 

from   the  analysis  of  all  runs  and  plots  we  have  the 
following  results: 

1.  Simulating  the  type  S#25#rr  and  S#5C#rr  model  (with 
r  =.25,  .50  and  .90)  it  was  possible  for  the  W  and  w"  to 
reach  the  steady  state  for  N=2000.  On  the  ether  hand  for 
the  type  S#25#rr  and  S#50#rr  (with  rr  =  .95,  .98)  it  was 
necessary  to  go  up  to  N=10000,  in  order  for  the  W  and  fl  to 
reach  the  steady  state.  Thus  we  see  that  the  high 
correlation  affects  the  choice  cf  the  value  of  N,  requiring 
N  to  increase  as  the  correlation  is  increasing.  But  not 
only  the  correlation  affects  the  choice  of  N.  The  traffic 
intensity  affects  it  luch  mere,  since  for  tne  type  S#95#rr 
and  S#99#rr  the  W  and  TJ  do  net  reach  the  steady  state  for 
N=10000  even  if  r  is  low.  Because  cf  the  CFQ-time 
requirements  we  restricted  ourselves  to  a  detailed  study  of 
the   particular   case   S#99#98   increasing   successively  the 

value  of  N  up  to  320,000,  where  the  W  and  W    appear   to   be 

close  to  the  steady  state.  Thus  we  may  ccnclude  that  high 
traffic  iutersity  and/cr  high  correlation  require  a  large 
value  of  N  in  order  to  achieve  steady  state.  Figures  23a 
through  23d,  24a-24c,  25  and  26a-26d  justify  the  abeve 
conclusicns.   In   figures  23a-23d,  where  we  are  dealing  with 


117 


th 

the  waiting   time  cf  the  2000   (w    )  arrival  for  the   case 

2000 


t=0.50,   r=0.25, 


we   can  see  that  »      has  converged  to  a 

2000 


value  of  about  6.3,  and  furthermore  a  straight  line   appears 

in    the    plot    under   EXPLT.    We   also   see   approximate 

convergence  for  the  correlated  gueue  for  the  H  ,    which  is 

2000 
presented   en   figures   28a-28c,  and  its  value  has  converged 

to  a  value  cf   about   3.03   (figure   23d)  .    Note   that   the 

introduction   of   correlation   into   the   service   time   has 

increased  the  average  waiting  time  from  2.5  tc   about   3.03, 

cr  by  approximately  2Q%. 

We  cannot  say  the  same  icr  figures  2Ua-24c  where, 
because  cf  the  high  ccrrelaticn  (r-0.98)  we  do  not  have 
convergence  cf  W  and  H,   although  N=10,000,  and   for    figure 

25  where,  because  cf  the  high  traffic  intensity  (t=0.S9)  we 
dc  not  have  convergence  of  W  and  ¥.  Figures  26a-26d  also 
show  us  the  ncn-ccnvergence  of  W,  H,  because  of  the  high 
correlation  and  the  high  traffic  intensity  (t=0.99, 
r=0.98) .  This  is  the  case  for  which  we  eventually  had  to 
carry  the  siaulaticn  out  to  about  N=500,000  to  observe  a 
stationary  gueue. 
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Figure    23a    -         QUEUE       WITH       EAR1    AUTOREGRESSIVE    SERVICE    TIME 
SEQUENCE    AND    POISSON    INPUT.        HISTOGRAM       OF    IHE    WAITING    TIMES 


H 


WITHOUT 


ZEROS 


FROM 


THE 


RUN 


20  00 

REPLICATIONS,    RX=.2,    RS=.  4.  No    OF    ZEROS  =  255. 


S#50#25;       m=500 
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U   U     H 


figure  23t  -  QUEUE  WITH  EAB1  AUTOREGBESSIVE  SERVICE  HflE 
SEQUENCE  ANE  POISSON  INPUT.  EXPCNENTIAI  PICT  (EXPIT)  CF  THE 
WAITING  TIMES  W      WITHOUT  ZEROS   EROS   TSE   RUN    S#50#25; 

2  00  0 

ffl=500    REPLICATIONS,    RX=.2,     RS=.4.  No    OF    ZERCS=255. 
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♦  CD 

♦  • 

«  O 


Figure  25c  - 
SEQUENCE   AN£ 


QUEUE  WITH  EAE1  A UTCREG EESSI V E   SERVICE   TIME 
PCISSON   INPUT.    NORMAL  PLOT  (NORMPL)  CF  THE 

FROM   THE   RUN    S#E0#25; 
No  OF  ZF5CS=255. 


WAITING  TIMES  W      WITHOUT  ZEROS 

200  0 

m=500  REPLICATIONS,  RX=.2,  RS=.4. 
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figure    24a    -         QUEUE    WIT3     EAR1    AUTOREG RESSIVE       SERVICE       TIME 
SEQUENCE       ANE    PCISSON    INPUT.       HISTOGRAM    OP    THE    WAITING    TIMES 


W  WITHOUT       ZEROS       FROM 

10000 

REPLICATIONS,    RX=2.5,     RS=10; 


THE         RUN  S#25#98; 

No    OF    ZEROS=354. 


m=500 
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Figure  24t  -  QUEUE  HUH 
SEQUENCE  AND  POISSON  INPUT 
WAITING   TIMES 


EAR1  AUTOREGRESSIVE  SERVICE  TIME 
EXPONENTIAL  PICT  (EXPLT)  CF  THE 


N        WITHOUT  ZEROS  FROM  TEE  RUN   S#25#98; 

lOOOO 

m=500  REPLICATIONS,  RX=2.5,  RS=10;   No  0?  ZEEQS=354. 
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Figure  24c  -  QUEUE  WITH  EAR1  A UTORSG RESSI VE  SERVICE  TIME 
SEQUENCE  ANE  POISSON  INPUT.  NORMAL  PLOT  (NORMPL)  CF  THE 
WAITING  TIMES  W       WITHOUT  ZEROS  FROM   THE   RUN    S#25#98; 

lOOOO 

m=500  REPLICATIONS,  RX=2.5,  RS=10;   No  OF  ZEE0S=354. 
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Figure    24e    -         QUEUE       WITH    EAR1     AUTOREGRESSI VE    SERVICE    TIMES 
AND    FOISSCN    INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND    AVERAGED 


WAITING 


TIMES 


¥ 


FROM 


THE 


100  00 

REPLICATIONS,    RX=2.5,     RS= 1 0 


RUN       S#25#98;       m=500 
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figure  24f  - 
AND    PCISSON 


QUEUE  WITH  EAE1  AUTOREGRESSI V E  SERVICE   TIMES 
INPUT.     EXPONENTIAL   ELOT   (EXPLT)    OP   THE 


CUMULATEE  AND  AVERAGED  WAITING  TIMES  H 

lOOOO 

S#25#98;     db=500    REPLICATIONS,    RX=2.5,    RS  =  10     . 
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Figure  2<4g  - 
AND   POISSON 


QUEUE  WITH  EAB1  A UTOREGRESSI V I    SERVICE   TIMES 
INPUT.    NORMAL  FLCT  (NOfiflPL)  CF  THE  CUMULATED 


AND  AVERAGED  WAITING  TIMES   TJ 

100  00 

m=5Q0  REPLICATIONS,  RX=2.5,  RS=10 


FROM   THE   RUN   S#25#98; 
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Figure    26a    - 
AND       P0ISS0N 


QUSOS    WITH     2AR1     AUTOREGRESSIVZ    SERVICE       TIMES 
INPUT.        HISTOGRAM       OF    THE    WAITING    TIMES     W 


10000 

WITHOUT    ZEROS    FROM    THE    S#99#98;     31  =  500    REPLICATIONS,     RX  =  2.97; 


RS  =  3 


Nc    OF    ZER0S=19. 
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Figure  26b  -  QUEUE  WITH  EAR1  A UTORSGRESSI VE  SERVICE  TIMES 
AND  POISSON  INPUT.  EXPONENTIAL  PLOT  (EXPLI)  OF  THE  WAITING 
TIMES  W       WITHOUT  ZEROS   FROM   THE    RUN   S#99#98;   m=500 

lOOOO 

REPLICATIONS,  RX=2.97;  RS=3  ;   No  OF  ZER0S=19. 
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Figure    26c    - 
AND       POISSON 


QUEUE    WITH     EAR1    A UTOREG RESSI VE    SERVICE       TIMES 
INPUT.  NORMAI       FIOT       (NORMPL)     OF    THE    WAITING 

FROM      THE         RUN       S#99#98;       m  =  500 


TIMES    W  WITHOUT    ZEROS 

10000 

REPLICATIONS,    RX=2.97;    RS=3     ;       No    OF    ZER0S=19. 
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2.   As   for  the  steady  state  of  the  value  cf  E[  W ]  and  E[w"] 
a   large   N   was   required  for  high  t  and/or  r,  also  for  the 

steady  state  of  the  distibution  cf  W  and  "J?   a   large    value 

of  N  was  required.  Informally  testing  their  distributions, 
ty  using  the  Plotting  Subroutines  and  the  values  of  their 
parameters  skewness  and  kurtosis,  we  may  conclude  that: 

(i)  .  The  distribution  of  W,  given  that  W>0,  has  an 
exponential  form,  if  the  plot  results  in  a  straight  line 
under  EXFLT  and  the  '/^and  JL parameters  have  the  values  of 
about  2  and  6  respectively  (see  figures  23a-23d) .  This  is 
true  only  if  the  steady  state  has  been  reached  while  it  is 
net  evident  hhat  their  distribution  is  if  the  steady  state 
has   not   been  reached.   Thus  observing  the  figures  24a-24d, 

where  the  W      has  not   reached  the  steady  state   yet,   we 
1C000 

can  obtain  nothing  about:  the  distribution  of  W. 


Note  that  for  the  M/M/ 1  queue  the  steady  state 
distribution  of  W,  given  that  W>0 ,  is  exponential.  The 
distribution  is  not  known  for  the  correlated  gueue  and  this 
was  one  cf  the  objectives  of  this  study. 

(ii)  .   The  distribution  of  W*   has  the  normal  form  if   we 

get    a   straight   line   under   NOBMPL   and  if  the  values  of 

skewness   and   kurtosis   are   around   fQ'.  The    ncrnality 

assumption   holds   for   the   steady   state  only,   while  the 

non-steady  state  distribution  looks  like  an  exponenticnal. 

Figures   27a-27c   give  us   the   feeling   that  w*      tas  an 
.  10000 

exponential  form,  since  a  straight  line  appears  under   EXPLT 

(see  figure  27b)  and  the  parameters  V^,  j)^  have  the  values 
2.3  and  7.3  respectively.  These  are  close  to  the  cctual 
values   2  and  6  of  the  exponential  distribution.   We  can  see 
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from  the  clefs  that  as  the   T?    gets   close   to   the   steady 

state,  it  leaves  the  exponentional  form  and  ccmes  closer  and 

closer  to  the  normal  form.    Figures  28a-28c,   23d   shew   us 

that  1?       has   already   reached  the   steady   state   and 
10000 

furthermore  its  distribution  has  eventually  taken  the  normal 

form,   since   a   straight  line  appears   under  NORMPL  (figure 

28c)  and  the  skewness  and  kurtosis  have  the  value  .7  and   .8 

respectively.      Compare    its   distribution    with   the 

distribution  of  w*      of  figure  27b  where  1  is  not   in   the 

10000 

steady   state.    Again  we  may  state  here  that  convergence  to 

the  steady  state  of  the  distribution  of  W   and  ¥    reguires 

much  larger  N  than  convergence  to  the  steady  state  of  their 
mean  values  for  the  same  t  and  r.  And  alsc  the  higher  t 
and/or   r   the   larger  N  should  be.  As  an  example  we  can  see 

that  even  thcugh  the   steady   state   (in   value)   of   TJ   for 

S#25#98   type  has  teen  reached  fcr  N=100Q0,  the  steady  state 

in  distribution  has   net   been   reached   yet.    See   figures 

2<4e-24g  where  we  can  otserve  that  the  distritution  of  15 

10000 
starts  to  leave  the  exponential   form   and   to   go   tc   the 

normal  form. 
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m  =  500 
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Figure    211    -         QUEUE       WITH    EAR1     AUTOREGRESSI VF,    SERVICE    TIMES 
AND       POISSON       INPUT.  EXPONENTIAL       PLCT        (FXPLI)        CF         THE 

CUMULATED       AND       AVERAGED       WAITING       TIMES 


"R  FROM    TEE    RUN 

2  000 


S#5C#9Q;     m=500    REPLICATIONS,    RX=.2;    RS=.4 
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Figure  27c  -  QUEUE  HIIH  EAR1  AUTOREGRESSI VE  SERVICE  TIMES 
AND  POISSCN  INPUT.  NORMAL  HOT  (NORMFL)  CI  THE  CUMULATED 
ANE  AVERAGE!  WAITING  TIMES  ¥      FROM  THE  RUN  S#50#90;  m  =  5GO 

2000 

REPLICATIONS,  RX=.2;  RS=. 4  . 
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Figure    28a    -         QUEUE    WITH    EAR1    AUTORSGRESSIVE    SERVICE       TIMES 
AND    POISSCN     INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND    AVERAGED 


WAITING       TIMES 


FROM 


20  00 

REPLICATIONS,    RX=.2;     RS=.4 


THE 


RUN 


S#50#25; 
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Figure    28t    -         QUEUE       KITH    EAR1     AUTOREGRESSI 11    SERVICE    TIMES 
AND       POISSON       INPUT.  EXPONENTIAL       PLCT        (FXPLT)        CF         THE 

CUMULATED       AN£       AVERAGED       WAITING       TIMES 


¥  FROM    THE    RUN 

2000 


S#5C#25;    m=5CC    REPLICATIONS,    RX=.2;    RS=.4 
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Figure  28c  -  QUEUE  KITH  EAR1  AUTOREGRESSI VE  SERVICE  TIMES 
AKD  PCISSCN  INPUT.  NORMAL  ELOT  (NORMPL)  01  THE  CUMULATED 
ANE  AVERAGE!  WAITING  TIMES  H  FROM  THE  RUN  S#50#25;  Dn  =  500 

2000 

REPLICATIONS,  RX=.2;  RS=.4  . 
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Continuing      the      analysis      of    the    results    we    present 


below    the   results    we      got      from      all      runs,      concerning 
values   of    E£"B]   for    various    values   of   t,    r,    N. 


the 


N 

2000 

2000 

2000 

10000 

1000G 


t 

25 

25 

25 

25 

25 


r  E[U]  (correlated)  E£¥]  (M/M/1) 

25  .3720  .3332 

50  .4326  .3332 

90  .7578  .3332 

95  1.0576  .3332 

98  1.7077  .3332 


2000 

2000 

2000 

10000 

10000 


.50 
.50 
.50 
.50 
.50 


25 

50 
90 

95 
98 


1.2106 
1.5811 
5.0785 
9.3196 
21.0542* 


.9985 
.9985 
.9985 
.9985 
.9985 


10000 
10000 
10000 
10000 

10000 


.95 
.95 
.95 
.95 
.95 


25 
50 
90 
95 
98 


23.5651 

33.3367 

123.9350* 

196.2460* 

342.9240* 


18.3528 
18.3528 
18.3528 
18.3528 
18.3528 


10000 
10000 
10000 
10000 
10000 


.99 
.99 
.99 
.99 
.99 


25 
50 
90 
95 
98 


62.3285* 

77.9370* 

193.9113* 

271.7713* 

426.6750* 


5  3.3371 
53.3371 
53.3371 
53.3371 
53.3371 


Note  than  these  values  have  been  scaled  from  the  values 
obtained  in  the  simulations  by  multiplying  by  RS.  Thus  for 
t  =  0.50  and  r=0.25  the  table  gives  2[  TJ  ]  (correlated)  as 
1.2106;  the  value  in  Figure  23d  is  3.026x0.4=1.2106  .  The 
values  marked  *  have  net  converged. 
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A  study  on  these   results  gives   us   the   following 
observations : 

1.  For  given  traffic  intensity  t  the  value  of  E[¥]  is  not 
constant,  as  it  is  for  M/M/1  queues,  but  it  depends  en  the 
value  of  the  correlation  parameter  r.  Thus  we  see  that  the 
value  of  E[w"  ]  increase,  as  the  value  of  r  becomes  larger 
and  larger.  Furthermore  from  figure  29a,  where  E[T»]  is 
plotted  versus  r  we  can  see  that  the  rate  cf  increment  is 
not  constant  (therefore  is  net  a  linear  function  of  r)  but 
it  increases  with  r  and  eventually  gees  up  infinitely  as  r 
goes  to  1. 

2.  For  given  correlation  r,  we  see  that  the  value  of  E[TJ] 
increases   nith   t,   the   sane  as  happens  in  the  M/M/1  gueue 

(see  figure  29b).  Eut  regardless  of  the  value  of  r,  E[ « ]  for 
M/M/1  gueue  is  always  less  than  for  the  EARMA  model.  Thus 
we  may  state  that  the  waiting  time  of  the  autocorr elated 
service  time  model  is  greater  than  the  corresponding  waiting 
time  of  the  M/M/1  model. 
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Figurs    29a    -         QUEUE    WITH    EAR1    AUTOREGSES5IVE    SEET7ICE 
HUES    AN£    POISSON    INPUT.       PLCT    OF    E[  ¥  ]    VERSUS    r. 
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figure    29b    -  QUEUE    WITH    EA&1    AUTOfiEGSESSIVE    SEFVICE 

IlflES    AN£    PCISSON    INPUT.       PICT    OF    E[ ¥ ]    VERSOS    t. 
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3.  The  values  of  E[  ¥"  ]  for  the  type  S#99#98  model 
increases  with  N  following  the  curve  in  figure  30  a  fact 
that  is  known  to  happen  in  the  M/M/1  queue  also.  We  can  see 
from  that  figure  that  the  rate  of  increment  comes  dcwn  as  N 
increases   and  hopefully  should  come  to  zerc   as   w*   reaches 

the   steady   state.   Comparing  then  the  M/M/1  queue  with  the 

correlated  queue  we  can  see   that  the  M/M/1  queue   converges 

faster   than  the  correlated  especially  when  r  is  high.   Thus 

we  can  see  from  figures  32a-32c  that   both   the   value    and 

distribution   of  the  1  for  M/M/1  queue  have  reached  the 

10000 

steady  state,  but  the  corresponding  T7      cf  the  correlated 

10000 
gueue    (figures    33a-d)     has    not    (at  least   in   its 

distribution).   (See  the  skewness  and  kurtcsis  values  for  W1 

EAR,   W2  EAR,  W3  BAR  and  W4  BAR  in  Figure  33d.   These  appear 

to  te  decreasing  tc  zero,  but   have   certainly   not   reached 

zero   by   N=10,000.)     Looking  at  figure  31c,  we  can  cbtain 

the  result  that  neither  U       fcr  the  correlated  queue,   or 

10000 

"R       fcr   the  unccrrelated  queue  has  reached  equilibrium. 
10000 

Ihis   is   because   the   skewness   and   kurtcsis   values   are 

definitely  net  close  tc  zero. 
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Figure    30    -         QUEUE    WITH    EAR1     AUTOREGRESSIVE    SERVICE    TIMES 
AND    POISSON    INPUT.       PlOr    OF    E[T7]    FROM    THE    RUN    S#99#98  ,    vs    N. 
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ISTOGRAU    OF    THE    CUMULATED    AND    AVERAGED 

FROM         THE         RUN       S#95#95;       m=500 
oo 
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Figure  31b  - 
AND   POISSON 


QUEUE   WITH  EAR1  AUTOREGRESSI V S  SERVICE  TIMES 
INPUT.     EXPONENTIAL   PLOT   (EXPLT)    OF    THE 


CUMULATED   AND   AVERAGED   WAITING   TIMES  ¥       FROM  THE  RUN 

10  0OO 

S#95#95;  m=500  REPLICATIONS,  RX=3.8;  RS=4  . 
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Figure  31c  -    QUEUE   WITH  EAR1  AUTOREGRESSIVE  SERVICE  TIMES 
AND  POISSON  INPUT.    NCRML  FLCT  (NORMPL)  OF   THE   CUMULATED 


AND   AVERAGE!   WAITING   TIMES 


1O00O 


FROM  THE  RUN  S#S5#95; 


m=500  REPLICATIONS,  RX=3.8;  RS=U 
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Figure    32a    -         M/M/1    QUEUE.        HISTOGRAM       OF    THE    CUMULATED    AND 
AVERAGED       WAITING       TIMES 


REPLICATIONS,     RX=4 ;     RS=8 


¥  FROM    THE    RUN    S#50#95;     m=500 

1000  0 
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Figure    32b    -         M/M/1       QUEUE.  EXPONENTIAL       ELOT     (EXPLT)        OF 

LND    AVERAGED    WAITING    TIM 
S#50#95;     m=500    REPLICATIONS,     RX=4;    RS=i 


THE    CUMULATED    ANC    AVERAGED    WAITING    TIMES    ¥  FROM    THE    RUN 

lOOOO 
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Figure  32c    -         M/M/1  QUEUE.    NORMAL  PLOT    (NCRMPL)    OF   THE 
CUMULATED   AND   AVERAGZD   WAITING   TIMES  w       FROM  THE  RUN 

lOOOO 

S#5C#95;  m=500  REPLICATIONS,  RX=4;  RS=3  . 


155 


•  •       IN 


• 

# 

r* 

• 

» 

» 

# 

1 
1 
1 

a 
m 

♦ 

•  * 

• 

•N 

•  • 

♦ 

•  • 

♦ 

1 
1 
1 

O 

1 

1 

r-( 

1 

* 

•  *  • 

•  • 

* 

in 

•  • 

■ 

-*j 

•  •  • 

W 

*  * 

♦ 

•  * 

♦ 

# 

» 

« 

• 

O 

* 

• 

?*• 

»  •  »  • 

rn 

*  »  ♦ 

• 

#  ♦  *■ 

* 

3> 

#«»#»•**** 


»*****.      o 


VSBsccssxflSvcccesscxxccssjcsc&CEsiicvcticrcciirerviisiii: 


ooooooo 
I  I   I 

vn  ^  *"  *"•  <■•">  ^  a? 


UUUJWU.UJ 


3  35  SE  :£  '-=  Z 

ai  i:  -c  -T.-X  «;  «*  tz 

W  3O0333ZD 

as  C003XIC 


CN(N*-0<"MfN 

oooooo 
III       ti 

a  laiifju'.r- 
»-o  T-*-  ^-3 


lilt       I 

r*ocf— u-\«— 

JT-g"t03-0 


n: 


U.J 


«0  MO* 


JJ  QU  CUUJlJ    U  JU 


^  qmuuj  -a 

LU  [■)XH-..J< 

c  *3 1- K  r;  ■  3  a: 


Figure    33a    -         QUEUE       WITH    EAR1     AUTOREGRESSIVE    SERVICE    TIMES 
AND    POISSOS    INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND    AVERAGED 


WAITING 


TIMES 


100  00 

REPLICATIONS,    RX=4 ;     RS=3 


FROM 


THE         RUN       S#50#95;       ie  =  5G0 
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Figure  33b  - 
AND   POISSON 


QUEUE   WITH  EAR1  AUTOREGRES SI VE  SERVICE  TIMES 
INPUT.     EXPONENTIAL.  PLCT   (EXPLT)    OF    THE 


CUMULATED   AMD   AVERAGED   WAITING   TIMES  ¥       FROM  THE  RUN 

10000 

S#50#95;  ni  =  500  REPLICATIONS,  RX  =  4;  RS=8  . 
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Figure  33c  -  QUEUE  WITH  EAR1  AUTOREGRESSI VE  SERVICE  TIMES 
AND  POISSON  INPUT.  NORMAL  ELOT  (NORMPL)  OF  THE  CUMULATED 
AND   AVERAGE!   WAITING   TIMES 


¥        FROM  THE  RUN  S#50#95; 

lOOOO 


01  =  500  REPLICATIONS,  RX=4;  RS=8 
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Another    result    obtained    from    the   simulation 

concerns  the  comparison  of  convergence  of  tt  and  w*.  Froi   the 

programs   we   have   run  we  observe  that  w"  converges  slightly 

faster  than  H  (at  least  for  low  t  and  r) .   Also  the  standard 

deviation  of  H    is  much  less  than  of  W.   In  figure  2Ud  we  can 

see  that  H  has  already  reached  the  steady  state  but  not 

10000 
*      .     Compare    also    the   values   of   their   standard 

10000 
deviations  (.005  versus  .2).   Figure  34  gives  the  results  of 

the   type  S#50#98  model  for  N=10000.    We  can  see  that  while 

1  is  already  close  to  the  steady  state,  the  fact  that  fc   has 

not   reached   it  yet  fcrces  us  to  continue  the  simulation  by 
increasing  the  value  of  N. 
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2.   £ueuj  with  Cross-Correlated  Service  Time 


The  model  for  the  cress-correlated  queue  is: 

<x.)  =  6. 

i        i 

{S    }         is  EAEMA(1,1)    over   {  E  ,    RX/RS  ,      RX/RS,  . 
i  i   i        i-1 


••  } 


That  is, 


I  #E.  +  A 


l     l 


where 


r  rA 


A   = 

i 


i-1 


w.p.  ^ 
w.p   (1-#) 

w.p.   r 


rA     +  g  RX/RS    w.p.   (1-r) 
i- 1     i 


Again  we  chose  the  values  .25,  .50,  .95,  .99,   for  t 

and   .25,   .50,  .90,  .95,  .98  for  r  and  we  run  a  total  cf  20 

runs  to  cover  all  cc urinations  of  t  and  r,  naming   each   run 
as  SC#tt#rr. 

Analyzing  the  results  obtained  from   the   simulation 
we  cbservs  the  following: 

1.   The    value    of  the   traffic   intensity   t   and   the 

correlation  r  is  the  main  factor  for  convergence  of  W  and  fl. 

That   is,   high    traffic   intensity  and/or  high  correlation 

require  N  to  be  large  in  order  for  W   and  H      to  reach   the 

steady   stats.    Thus   simulating   the  models  with  t<.50  and 

r<.95  for  N=10000  it  was  possible  for  *  and  H   to  reach   the 

steady   stats.     In  figures  35a-35d,  36a-36c  we  can  sss  the 

justification  of  that   conclusion.    The   values  .1076   and 

.0957   of  ~A  ,   and  ¥      represent  the  steady  state  of 

10000         10000 
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the  SC#5G#25  model.  On  the  ether  hand,  regardless  of  the 
value  of  t  the  models  with  r>.98  did  nor  reach  the  steady 
state  for  N=10000  and  regardless  of  the  value  of  r  the 
ircdels  with  t>.95  did  net  reach  the  steady  state  for 
N  =  1 G  #  00  0  .  figure  37  gives  us  the  non-convergence  of  the 
model  SC#99tf25.  Because  of  the  high  traffic  intensity  and 
in  order  to  get  the  steady  state  of  that  model  we  have  to 
continue  simulating  by  increasing  the  value  cf  N. 
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Figure    35a    -         QUEUE         WITH 
SEQUENCE    AND    POISSON    INPUT. 


N 


WITHOUT 


ZERCS 


10000 

REPLICATIONS,    RX=2;    RS=4 


CROSS-CORRELATEE       SERVICE  TIME 

HISTOGRAM    OF    THE    WAITING  TIMES 

FROM         THE      RUN       SC#50#25;  m=500 

No    OF    ZEROS=279. 
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r*      ?•    ^ 


Figure  35t  -  QUEUE  WC  TH 
SEQUENCE  AND  POISSON  INPUT. 
WAITING   TIMES 


CROSS-CORRELATEE   SERVICE   TIME 
EXPONENTIAL  PLCT  (EXPLT)  CE  THE 


W       WITHOUT  ZEROS  FROM  THE  RUN   SC#50#25; 

10000 

m=500  REPLICATIONS,  RX=2;  RS=4  .    No  OF  ZERCS=279. 


165 


*  • 


•  <7> 

♦  • 


#  3 


Figure  35c  -  QUEUE  WITH 
SEQUENCE  AND  POISSON  INPUT. 
WAITING   TIMES 


CROSS-CORRELATEE   SERVICE   TIME 
NORMAL   PLOT    (NORMPL)    OF   THE 


W       WITHOUT  ZEROS  FROM  THE  RUN   SC#50#25; 

lOOOO 

ffl=500    REPLICATIONS,    RX=2;     RS=4     .  No    OF    ZERCS=279. 
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2.  High  trafic  intensity  and  correlation  affects  the 
convergence  of  the  distribution  of  W  and  ¥  also. 
Furthermore  the  steady  state  of  their  distributions, 
reguires  larger  N  than  is  reguired  for  the  steady  state  of 
their  expected  values.  In  examining  that  conclusion,  look 
for  example  at  the  results  cf  the  simulated  model  SC#50#98, 
where   although    W   and  W"   have  already  reached  the  steady 

state  for  N=1Q000,   their   distribution   has   not   converged 

yet.     The   above   results  have  been  obtained  from   Figures 

58a-38c,  5Sa-39d  wnere  W       and  ¥      are  presented.    We 

10000       10000 

can  see  from  tnese  figures  (38a-38c)  that  although  the  value 

of  W      has  reached  the  steady   state   (see   Figure   39d) , 

10000 
its   distribution   cannot   be    characterized  since  the  plot 

under  EX£1T  (figure  38b)  is  not   a   straight   line   and   the 

parameters    ]jL=3.  34    and  ^  =  11.0   are   far   away   from  the 

desired  values  2  and  6  .    I  he  same  results  can  be   obtained 

from  the  figures  39a-3Sc  which  deal  with  ¥ 

10C00 
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Figure    36a    -         QUEUE  WITH       CSOSS-CORRELATEE       SERVICE       TIME 

SEQUENCE    AND    POISSON    INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND 

AVERAGED       WAITING      TIMES    ¥  FROM    THE    RUN    SC#50#25;     m=500 

10  000 

BEPLICATICNS,    RX=2 ;    RS=4    . 
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Figure    36t    -        QUEUE         WITH 
SEQUENCE    AND    P0ISS0N    INPUT. 


CROSS-CORRELATEE       SERVICE       TIME 
EXPONENTIAL       PICT        (EXPLT)        OF 


THE    CUMULATED    AND    AVERAGED    WAITING    TIMES    H  FROM    THE    RUN 

10  000 

SC#50#25;    m=500    REPLICATIONS,    RX=2;    RS=4    . 
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Figure    36c    -         QUEUE         WITH 
SEQUENCE    AND    POISSON    INPUT. 


CROSS-CORRELATEE       SERVICE       TIME 
NORMAL    PLOT        (NGRMPL)        OF      THE 


CUMULATED       AND       AVERAGED       WAITING      TIMES    ¥  FROM    TEE    RUN 

lOOOO 

SC#50#25;    m=500    REPLICATIONS,    RX=2;     RS=4     . 
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Figure  38a  -  QUEUE  WITH  CROSS-CORRELATED  SERVICE  TIME 
SEQUENCE  ANC  POISSON  INPUT.  HISTOGRAM  Or  THE  WAITING  TIMES 
W  WITHOUT       ZEROS      FROM      THE         RUN  SC#50#98;  m=500 

l  0000 

REPLICATIONS,     RX=2;     RS=4     .  No    OF    ZEROS=245. 
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Figure  38b  -  QUEUE  WIT H  CROSS-COBR ELATED  SERVICE  TIME 
SEQUENCE  AND  POISSON  INPUT.  EXPONENTIAL  PLOT  (EXPLT)  CE  THE 
WAITING  TIMES  W       WITHOUT  ZEROS  FROM  THE   RUN    SC#50#98; 

lOOOO 

m=500  REPLICATIONS,  RX=2;  R5=4  .    No  OF  ZERCS=2U5. 
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Figure  38c  -  QUEUE  WIT H  CfiOSS-CORRELATEE  SERVICE  TIME 
SEQUENCE  ANC  POISSON  INPUT.  NORMAL  PLCT  (NORBPL)  CF  THE 
WAITING  TIMES  W       WITHOUT  ZEROS  FROM  THE   RUN    SC#50#98; 

1  00  00 

m  =  5  00  REPLICATIONS  ,  RX  =  2;  RS  =  U  .    No  OF  ZEECS  =  245. 
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Figure    39a    -         QUEUE       WITH       CROSS-CORRELATED       SERVICE         TIME 
SEQUENCE    AND    P0ISS0N    INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND 


AVERAGED    WAITING       TIMES 
REPLICATIONS,    RX=2 ;     RS=4 


1  0  000 


FROM       THE       SC#50#98;       m=500 
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Figure  39b  -    QUEUE   WITH   CROSS-CORRELATED   SERVICE    TIME 
SEQUENCE   ANC   POISSON   INPUT.    EXPONENTIAL  ELOT  (EXPIT)  OF 


THE  CUMULATEE  AND  AVERAGED  WAITING   TIMES 
SC#50#98;  a=500  REPLICATIONS,  RX=2;  RS=4 


lOOOO 


FBOE   THE 
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figure    39c    -         QUEUE       WITH       CROSS-CORRELATED       SERVICE  TIME 

SEQUENCE       ANE       POISSON       INPUT.  NORMAL    PLOT     (NORMPL)     CE    THE 


CUMULATED       AND       AVERAGED       WAITING       TIMES 
SC#50#98;    a=500    REPLICATIONS,    RX=2;     RS=4 


lOOOO 


FROM 


THE 
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3.   To  tast  the  distribution  of  H  and  U    we  also   usee   the 

plotting  subroutines  and  the  values  of  their  skewness  and 
kurtosis  parameters  as  we  did  with  the  gueue  with  dependent 
service. 

The   results   obtained   from   the   analysis   are   as 
fellows: 

(i) .  Given  that  W>0,  we  may  say  that  H  has  an 
exponential  form,  since  the  plots  under  EXPLT  Subroutine 
result  in  a  straight  line  and  the  values  of  /v.  and  X^are 
around  their  actual  values  2  and  6  respectively.  For  example 
the  SC#50#25  model  (which  is  represented  by  the  figures 
35a-35d)  ,  where  the  steady  state  has  been  reached,  gave  us  a 
straight  line  and  the  values  1.8  and  4.7 
for  ^L  and  ^respectively.  We  have  to  state  here  that  the 
exponential  distribution  comes  only  if  W  has  reached  the 
steady  state,  and  we  cannot  have  the  same  result  if  the  N  is 
still  in  a  transient  state. 


(ii)  For  the  distribution  of  17  we  can  say  that  the  \ 
normality  assumption  is  an  appropriate  one,  since  the 
graphs  obtained  under  NORMEL  result  in  a  straight  lire  and 
the  value  of  skewness  and  kurtosis  is  around  '0'.  The 
results,  for  example,  obtained  from  the  SC#25#25  model  gave 
us  a  straight  line  under  NOEMEL  and  the  values  .17  and  -.05 
for  the  skewness  and  kurtosis  respectively  (see  figures 
40a-40c)  .  The  normality  assumption  holds  if  the  S  has 
reached  the  steady  state.  At  the  begining  of  the  transient 
state  it  is  net  possible  to  characterize   its   distribution, 

while   when   TJ   is   in  a  transient  state  we  can  see  that  its 

2 

distribution  takes  an  exponential  or  a   X    fcrm.    This   is 

left   little   by  little  as  TJ  comes  close  to  the  steady  state 
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and  then  eventually  it  takes   the   normal   fcrm.   The   model 
SC#25#95   is   a   representative   example  of  that  result.   In 

that  model,  although  w*  has  net  reached  the  steady  state   yet 

(N=  10000 )  we  can  observe  the  values  of  skewness  having  as 
3.5,  1.8,  1.6  1.5  and  the  values  of  kurtosis  as  24.7,  7.2, 
5.7,  4.3  foi  N=2500,  5000,  750C,  10000  r espectivelly .  lhat 
is  the  values  start  frcm  a  high  level,  then  ccme  down,  they 
pass  frcm  the  actual  values  of  ~i  ,  ^L  fcr  an  exponential 
model  and  then  continue  decreasing  and  hopefully  in  the 
steady  state  they  should  reach  a  value  arcurd  '0'  where  the 
value  of  skeviness  and  kurtosis  fcr  the  normal  distribution. 
The  same  assumptions  are  obtained  analyzing  the  plcts  of 
that  model  also.  That  is  for  N=2500  the  plots  are  not 
linear  under  both  EXPLT  and  N0RHPL,  for  N=5000  a  straight 
line  appears  under  EXPLT,  for  N=7500  and  10000  the  plot 
leaves  gradually  the  straight  line  under  EXEL1  and  begins  to 
give  us  a  straight  line  under  NOBMPL. 
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Figure    40a    -         QUEUE       Wir H       CROSS-CORSELATED       SERVICE         TIME 
SEQUENCE    AND    P0IS30N    INPUT.  HISTOGRAM    OF    THE    CUMULATED    AND 


AVERAGED    WAITING       TIMES      ¥ 


10  000 


FROM      THE       SC#25#25;       m=5G0 


REPLICATIONS,    RX=1.5;     HS=6 
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Figure    40t    -         QUEUE       WITH       CROSS-CORRELATED       SERVICE  TIME 

SEQUENCE       AND       POISSON       INPUI.  EXPONENTIAL    PLOT     (EXPIT)     OF 


THE    CUMULATED    AND    AVERAGED    WAITING       TIMES       U 


SC#25#25;    fn  =  5 0 0    REPLICATIONS,    RX  =  1.5;    RS=6 
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Figure   40c    - 
SEQUENCE       AND 
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4.   The  following  values  for  W"  have  been  attained  frcm  the 
runs  with  various  values  of  t  and  r  and  for  N=10000 


N 
1000C 
10000 
10000 
10000 
10000 

10000 
10000 
10000 
10000 
1000C 

10000 
10000 
1000C 
10000 
1000C 

10000 
10000 
10000 
10000 
10000 


t 

.25 
.25 
.25 
.25 

.25 

.50 
.50 
.50 
.50 
.50 

.95 
.95 
.95 
.95 

.95 

.99 
.99 
.99 
.99 
.99 


r 

.25 
.50 
.90 
.95 
.98 

.25 
.50 
.90 
.95 
.98 

.25 
.50 
.90 
.95 
.98 

.25 
.50 
.90 
.95 

.98 


B£¥] 

E[*]<iViV1) 

.1169 

.3332 

.  1724 

.3332 

.3538 

.3332 

.4075 

.3332 

.4679 

.3332 

.3831 

.9985 

.5648 

.9985 

1.4589 

.9985 

2.  1109 

.9985 

3.7196* 

.9985 

10.4672 

18.3528 

13.4784* 

16.3528 

29. 1312* 

16.3528 

75. 1180* 

16.3529 

14  5.3484* 

18.3528 

35.5986* 

53.3371 

43. 1616* 

53.3371 

98.9556* 

53.3371 

139.0167* 

53.3371 

219.7920* 

53.3371 

Analyzing  the  above  results  we  may  state  conclusions 
as  follows: 

(i)  .  Fcr  a  given  traffic  intensity  t  the  value  of  TJ  is 
net  constant,  but  depends  on  the  value  of  ccrrelation  r ,  as 
happens  in  the  queue  with  dependent  service.  We  observe 
also  that  "3  increases  as  r  increases,  but  the  rate  of 
increment  is  larger  as  r  goes  to  1,  and  very  small   for   low 
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values  cf  r;  therefore  ¥  is  not  a  linear  function  cf  r. 
Figure  41  gives  the  plots  of  ¥  versus  r  for  t  =  .25,  .50, 
.95,  .98  and  also  for  the  M/M/1  queue  (r=0)  .  Furthermore  we 
can  observe  that  the  value  of  ¥  for  the  correlated  gueue, 
comparing  it  with  the  result  for  the  M/M/1  queue,  is  not 
always  greater  as  in  the  case  cf  dependent  service,  but 
depends  en  the  value  of  r.  It  is  less  for  r<.50  and  greater 
for  r>.90 

(ii)  .  Fcr  a  given  correlation  the  value  of  ¥  increases 
as  t  increases.  Furthermore  high  values  of  t  cause  high  rate 
of  increment,  which  becomes  larger  as  the  correlation 
increases.  figure  42  gives  the  plots  cf  ¥  versus  t  for 
r  =  0,  .25,  .50,  .90,  .95,  and  .98  .  Observing  that 
figure  we  can  see  also  that  fcr  a  given  t  the  order  of 
increment  of  ¥  is  for  r  =  .25,  .50,  C,  .90,  .95,  .98 
respectively.  Thus  we  can  conclude  that  the  waiting  time 
fcr  the  cross-correlated  queue  is  generally  greater  than  the 
waiting  time  for  M/M/1  queue  when  the  correlation  is  high. 

5.  The  last  result  obtained  from  the  simulation  of  that 
model  was  the  convergence  of  W  and  ¥.  (Ihese  statistics 
alsc  converge  in  value  and  in  distribution  slower  than  the 
same  statistics  in  M/M/1  queue).  The  mcdel  SC#25#98  reveals 
that  ¥  converges  faster  than  W  (see  figure  43) .  Eut  the 
model  SC»99#50  gave  us  the  conclusion  that  H  converges 
faster  than  1  (see  figure  44) .  Thus  we  can  state  here  that 
fcr  low  values  of  t  ¥  converges  faster  than  H  and  fcr  high 
values  of  t  the  opposite  is  true. 
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Figure    41    -         QUEUE    WITH    CROSS-CORRELATED    SERVICE    TIMES    AND 
ECISSON    INPUT.       PLOT    OE    E[tf]    VERSUS    r. 
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Figure    42_  -         QUEUE    WITH    CROSS-CORRELATED    SERVICE    TIMES    AND 
PCISSON    INPUT.        PLOT    OF   E|]¥]    VERSUS    t. 
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